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ABSTRACT 


Upper  and  lower  bounds  were  determined  for  a  variation  of  Schmidt's  statistic 
using  ImhofF s  distribution  for  quadratic  forms  in  normal  variables.  This  statistic  is 
able  to  detect  a  fourth  order  autoregressive  disturbance  of  the  form:  Ct  =  PjCj-.j  + 
P4£t-4  +  Ht  general  linear  model  Y  =  XP  +  c. 

To  correct  for  this  disturbance  and  thus  yield  efficient  regression  estimates,  a 
data  transformation  was  derived  using  the  inverse  of  the  variance-covariance  matrix  as 
defined  by  Siddiqui. 
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I.  INTRODUCTION 


A.  BACKGROUND 

Extensive  work  has  been  accomplished  recently  in  the  area  of  modeling  and 
predicting  quarterly  overhead  costs  for  aircraft  manufacturers.  Overhead  costs  are 
generally  predicted  utilizing  estimated  overhead  rates  applied  to  labor  hours  or  costs 
over  functional  categories  such  as  manufacturing  or  engineering.  Changes  in  operating 
rates  cause  overhead  rate  changes  which  may  be  observed  only  after  a  significant  lag  in 
time.  Recent  work  done  in  this  area  has  been  to  estimate  total  overhead  as  a  function 
of  the  number  of  direct  manufacturing  personnel  [Ref.  1]. 

Data  collected  for  purposes  of  estimating  overhead  costs,  since  it  is  of  a  time 
series  nature,  can  be  expected  to  exhibit  some  degree  of  autocorrelation.  Specifically, 
data  collected  for  recent  purposes  has  exhibited  first  order  autoregressive  AR(1),  fourth 
order  autoregressive  AR(4),  or  a  combination  of  these  processes.  An  AR(1)  process 
occurs  when  the  errors  in  adjacent  time  periods  are  related.  This  type  of  relationship 
would  be  expected  in  yearly  cost  and  operating  data.  A  special  form  of  the  fourth 
order  autoregressive  would  be  expected  in  data  of  a  seasonal  nature,  i.e.  quarterly 
observations.  Errors  in  this  special  case  of  the  AR(4)  process  would  be  related  to 
errors  which  occured  four  quarters  previously. 

The  method  of  ordinary  least  squares  (OLS)  is  used  to  regress  the  data  and 
estimate  overhead  costs.  Ordinary  least  squares'  validity  is  based  on  certain  key 
assumptions: 

•  Errors  are  distributed  independently  of  the  explanatory'  variables  with  zero  mean 
and  constant  variance. 

•  Successive  errors  are  independent  of  each  other. 

With  autocorrelation  present,  assumption  two  is  violated  and  the  ordinary  least 
squares  procedure  breaks  down  at  three  points: 

•  Estimates  of  regression  coefficients,  though  unbiased  are  not  efficient. 

•  The  usual  formula  for  the  variance  of  an  estimate  no  longer  applies  and  is  liable 
to  seriously  underestimate  the  true  variance. 

•  t  and  F  distributions,  used  for  making  confidence  statements  are  no  longer  valid 
(Ref.  2). 

Since  the  OLS  procedure  breaks-down  with  the  presence  of  autocorrelation,  it  is 
important  that  its  presence  be  detected  and  the  disturbance  in  the  data  subsequently 
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corrected.  Without  accounting  for  these  disturbances  in  the  data,  estimates  of  the 
regression  coefficients  will  not  have  minimum  variance. 

The  Durbin-Watson  statistic  is  utilized  to  detect  the  presence  of  AR(1)  behavior. 
Additionally,  it  has  been  generalized  to  detect  higher  order  autoregressive  processes. 
Durbin  and  Watson's  work  was  based  on  the  findings  of  T.W.  Anderson  which  showed 
that  the  statistic 

e'Ae  /  e'e,  (eqn  1.1) 

where  e  is  the  column  vector  of  residuals  from  a  regression  and  A  is  a  certain  real 
symmetric  matrix,  provides  a  test  that  is  uniformly  most  powerful  against  certain 
alternative  hypotheses  [Ref.  2). 

The  general  model  upon  which  Durbin  and  Watson's  work  was  based  is 


(eqn  1.2) 


where 


st  =  etct_i  +  nt<  t=  i, . t 


(eqn  1.3) 


where  Xt  is  a  (1  x  K)  non-stochastic  matrix,  P  is  a  (K  x  1)  vector  and  i  is  one  for  an 
AR(1)  process  and  four  for  an  AR(4)  process  or  a  mixture  of  both  processes.  The 
error  component,  £t,  is  specified  by  the  form  above  and  r|t  is  distributed  with  zero 
mean  and  constant  variance. 

Typically  an  ordinary  least  squares  regression  is  performed  on  selected 
independent  variables.  If  the  presence  of  autocorrelation  is  suspected,  the  data  may  be 
tested  for  a  particular  AR  process  utilizing  the  appropriate  form  of  the  Durbin-Watson 
statistic.  For  an  AR(1)  process  the  form  of  the  Durbin-Watson  statistic  is 


(eqn  1.4) 


where 


(eqn  1.5) 


and 
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(eqn  1.6) 


The  estimator  of  P  from  the  ordinary  least  squares  regression  is  $  The  null  hypothesis 
for  the  test  is  that  of  zero  autocorrelation  in  the  residuals  against  an  alternative 
hypothesis  that  a  first  order  autoregressive  process  exists. 

If  the  presence  of  an  autoregressive  process  is  detected  in  the  regression  residuals, 
the  linear  model  must  be  reestimated  after  transforming  the  original  data.  For  an 
AR(1)  process,  the  independent  variables  are  transformed  as 


X,"  -  Xt(l  -  p,2)1/2  t-1 


(eqn  1.7) 


and 

Xt  =  Xj  -  t=  2,...T.  (eqn  1.8) 

In  a  similar  fashion,  the  dependent  variables  are  transformed  as 

Yt*  =  Yt(l  -Pl2)1/2  t=l  (eqn  1.9) 

and 

Yt*  =  Yt  -  PjYt_!  t  =  2,...T.  (eqn  1.10) 

In  order  to  accomplish  this  transformation  an  estimate  for  pj  is  required. 
Though  this  parameter  can  be  estimated  in  several  ways,  the  most  popular  method  due 
to  simplicity  of  calculation  is 

Pj  =  1  -  .5dj,  (eqn  1. 1 1) 

where  dj  is  given  in  equation  1.4  .  In  addition  to  the  ease  of  calculation,  this  estimator 
performs  as  well  as  more  complex  estimators  which  are  available  |Ref.  3]. 

After  the  data  is  transformed,  generalized  least  squares(GLS)  is  used  to 
reestimate  model  parameters.  After  GLS  is  performed,  the  model  is  rechecked  for  the 
presence  of  autocorrelation  using  the  Durbin-Watson  statistic. 
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B.  FOCUS  AND  SCOPE  OF  RESEARCH 

As  previously  mentioned,  regression  data  is  at  times  influenced  by  multiple 
autoregressive  processes.  In  overhead  cost  estimating,  both  an  AR(1)  and  AR(4) 
process  have  been  observed  to  simultaneously  influence  certain  data  sets  [Ref.  1],  In 
the  linear  least  squares  model 


Yt  =  Xtp  +  £t,  (eqn  1.12) 

suppose  the  errors  are  serially  correlated  and  of  the  form 

£t  =  0l£t-l  +  64£t-4  +  nt>  (ecln  113) 

where  the  effects  of  the  second  and  third  prior  quarters  are  negligible  compared  to  the 

effect  of  the  most  recent  and  year-earlier  quarters.  The  rjt's  in  the  model,  as  in  the 
previous  version,  are  independent  and  distributed  with  zero  mean  and  constant 
variance. 

It  is  advised  [Ref.  4:  p.  211]  that  if  this  special  form  of  the  AR(4)  process  exists, 
that  it  be  tested  for  by  a  two-step  procedure.  In  this  procedure,  after  initially 
performing  a  OLS  regression,  the  residuals  are  tested  for  an  AR(1)  disturbance  utilizing 
the  Durbin-Watson  statistic  for  a  first  order  process.  If  the  influence  of  an  AR(1) 
process  is  detected,  a  value  for  Pj  is  estimated  and  the  original  data  is  transformed. 
GLS  is  then  performed  on  the  transformed  data  and  the  residuals  are  tested  for  the 
influence  of  a  fourth  order  autoregressive  process  by  means  of  Wallis'  test  [Ref.  5].  If 
an  AR(4)  process  is  present,  a  value  for  p^  is  estimated  and  the  data  is  transformed  a 
second  time.  At  this  point,  the  procedure  is  repeated  to  determine  if  either  form  of 
autocorrelation  is  still  present. 

Unless  the  individual  performing  the  test  has  a  priori  knowledge  that  both  an 
AR(1)  and  AR(4)  process  exist  within  the  data,  it  is  very  likely  that  one  or  the  other 
would  be  overlooked.  Additionally,  peculiarities  in  the  data  may  preclude  the  detection 
of  one  or  the  other  processes.  If  this  occurred,  the  GLS  regression  estimates  as 
previously  discussed  would  be  inefFicient.  In  order  to  avoid  the  problem  above,  it 
would  be  desirable  to  have  a  procedure  to  detect  simultaneously  both  an  AR(1)  and 
AR(4)  process  if  they  existed  in  the  data.  This  process,  which  is  a  special  case  of  an 
AR(4)  process  shall  be  specifically  referred  to  as  an  AR(1,4)  process.  In  addition,  if 
the  AR(1,4)  process  exists  then  there  must  be  a  way  to  transform  the  data  and 
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calculate  estimates  for  the  parameters  of  the  process.  This  proposed  procedure  for  the 
AR(1,4)  process  would  prevent  oversights  of  existing  processes  and  would  save  time 
expended  in  having  to  correct  for  the  AR(1)  process  and  then  the  AR(4)  process. 
Most  importantly  though  it  would  insure  that  if  the  AR(1,4)  process  was  present  in  the 
data,  it  would  be  corrected.  In  this  way  any  peculiarities  in  the  data  which  might 
prevent  detection  and  correction  of  either  the  AR(1)  or  AR(4)  process  would  be 
avoided. 

As  previously  discussed,  with  an  autoregressive  disturbance  present  in  the  data, 
the  OLS  assumption  of  independent  error  terms  is  violated.  This  is  related  to  problems 
in  the  error  covariance  matrix  which  is  denoted  as 

E(ee')  =  (eqn  1.14) 

where  e  is  a  vector  of  regression  residuals.  The  problem  in  the  error  covariance  matrix 
is  that  there  exist  off  diagonal  terms  not  equal  to  zero.  This  condition  can  indicate 
that  some  form  of  autocorrelation  exists.  An  additional  problem  in  the  error 
covariance  matrix  is  that  the  diagonal  elements  are  not  equal.  This  condition  indicates 
that  the  error  terms  are  not  distributed  with  a  constant  variance.  Errors  distributed  in 
this  manner  cause  estimates  of  regression  coefficients  to  be  inefficient  though  unbiased. 
If  the  error  covariance  matrix  equals 

(eqn  1.15) 


instead  of 


(eqn  1.16) 


A 


where  I  is  an  identity  matrix,  the  regression  coefficient  P  is  correctly  defined  by 


(eqn  1.17) 


(X'vj/^XrX'y^Y. 


A 


The  general  procedure  [Ref.  6:  p.  440]  used  to  attain  an  efficient  estimate  for  P  is: 

•  Find  a  matrix  P  such  that  P'P  = 

•  Using  the  matrix  P,  transform  the  original  data  set  where 
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(eqn  1.18) 


Y*  =  PY 
and 

X*  =  PX.  (eqn  1.19) 


•  Perform  GLS  on  the  transformed  model 

Y*  =  X*p  +  e*  (eqn  1.20) 

where 

e*  =  Pe  (eqn  1.21) 

to  estimate  the  regression  coefficients  where 
/\  *  *  1  * 

p  =  (X  'X  )-‘X  Y  .  (eqn  1.22) 


C.  RESEARCH  QUESTIONS 

Following  the  procedure  as  outlined  above,  this  paper  will  address  the  AR(1,4) 
process.  Specifically,  this  thesis  will  attempt  to  develop: 

•  A  statistic  capable  of  detecting  the  AR(1,4)  disturbance. 

•  A  P  matrix  for  a  transformation  to  account  for  the  AR(1,4)  disturbance. 

•  A  method  to  estimate  parameters  required  for  the  transformation. 

D.  RESEARCH  METHODOLOGY 

In  reaching  the  first  objective,  the  statistic  will  be  calculated  by  estimating  the 
distribution  of  the  ratio  of  two  quadratic  forms  in  normal  variables.  This  method  is 
based  on  ImhofF s  [Ref.  7]  technique  and  is  outlined  by  Koerts  and  Abrahamse  [Ref.  8]. 
These  values  shall  be  calculated  for  an  extension  of  Schmidt's  statistic  [Ref.  9]  to  the 
AR(1,4)  process.  This  statistic  was  developed  in  a  manner  similar  to  Durbin  and 
Watson's  and  is  essentially  based  on  their  work.  The  power  of  this  statistic  to  detect 
the  AR(1,4)  process  will  be  compared  to  the  current  sequential  method  of  testing. 
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To  account  for  the  autoregressive  disturbance,  it  is  required  that  there  be  a  P 
matrix  capable  of  transforming  the  data  with  the  property  that 


P'P  =  y’1. 


(eqn  1.23) 


The  inverse  of  the  y  matrix  for  the  fourth  order  process  was  derived  by  Siddiqui 
[Ref.  10].  With  the  form  of  the  y  inverse  matrix  known,  a  P  matrix  with  the  above 
stated  property  can  be  found  using  a  method  proposed  by  Beach  and  MacKinnon 
[Ref.  11].  Both  the  y  inverse  matrix  and  P  matrix  are  expressed  in  terms  of  0's  or 
autoregressive  parameters  from  the  transfer  function 


(eqn  1.24) 


Utilizing  the  Yule-Walker  equations  [Ref.  12:  p.  55],  the  autoregressive  parameters 
from  this  transfer  function  were  expressed  in  terms  of  the  autocorrelation  coefficients 
or  p's.  Once  the  0's  were  determined  in  terms  of  the  p's,  substitutions  were  made  into 
the  derived  P  matrix.  The  resulting  P  matrix  can  be  used  to  transform  an  AR(4) 
process  as  represented  by  the  error  term  above.  Since  this  error  term  is  not 
representative  of  the  AR(1,4)  model,  corrections  were  made  to  the  P  matrix  by  setting 
appropriate  values  of  p  equal  to  zero.  For  Schmidt's  statistic,  which  will  be  introduced 
in  Chapter  II,  p->  and  p-j  were  set  equal  to  zero.  The  resulting  form  of  the  P  matrix  is 
that  required  to  transform  data  with  an  AR(1,4)  disturbance. 

With  the  P  matrix  defined,  the  next  task  is  to  find  estimates  for  pj  and  p^.  Since 
the  method  outlined  above  is  generalized  least  squares,  one  approach  for  estimating  the 
values  of  pj  and  p^  is  by  the  sample  correlation  coefficients 


rs  -  I  etet-s/I  et2  s=1>4 


(eqn  1.25) 


where  et  is  a  vector  of  least  squares  residuals.  Since  estimators  derived  in  this  manner 
are  subject  to  considerable  small-sample  bias,  the  approach  which  will  be  developed 
will  be  to  estimate  pj  and  p^  by  direct  utilization  of  the  test  statistic.  This  method  is 
analogous  to  Durbin  and  Watson's  method  where 


Pj  =  1  -  .5dj 


(eqn  1.26) 
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and  dj  is  the  Durbin-VVatson  statistic  for  an  AR(1)  process. 

After  developing  the  above  outlined  procedure  and  generating  the  bounds  for 
Schmidt's  statistic  a  major  question  which  must  be  answered  is:  How  well  do  the 
procedures  work  in  comparison  to  each  other?  In  order  to  make  the  comparison, 
simulations  will  be  used  which  utilize  least  squares  residuals  generated  for  various 
values  of  p  j  and  p_|  using  a  formula  for  the  error  of 

ct  =  Plct-1  +  P4et-4  +  nt  (ecln  *-27) 

where  T|t  is  distributed  normal  with  zero  mean  and  specified  variance  T  .  With 
residuals  calculated  as  such,  the  Schmidt,  Durbin-Watson  and  Wallis  statistics  will  be 
calculated.  In  addition  to  the  Durbin-Watson  statistic,  after  appropriately 
transforming  the  residuals,  the  Wallis  statistic  will  be  used  in  the  two-step  method  to 
test  for  the  presence  of  an  AR(4)  disturbance.  For  various  sets  of  residuals  generated 
with  different  variances  for  the  normal  error  term,  the  power  of  the  statistics  to  detect 
the  AR(1,4)  disturbance  will  be  determined.  Additionally,  the  number  of  observations, 
T,  will  be  varied  by  varying  the  number  of  residual  terms.  Thus  the  power  of  the 
various  statistics  to  detect  the  AR(1,4)  process  will  be  determined  for  different  values 
of  T,  pj  and  p_j  and  x  .  Residuals  used  in  the  above  calculations  will  also  be  used  to 
calculate  estimates  of  pj  and  p^.  Utilizing  mean  square  error,  calculated  from 
simulation  data,  the  efficiency  of  derived  formulas  to  estimate  Pj  and  p^  will  be 
determined.  Specifically,  the  ability  of  both  procedures  to  detect  the  AR(1,4) 
disturbance  will  be  determined  for  a  best  case  with  T=  100  and  x  —  .01  and  a  worst 
case  where  T  =  20  and  T2  =  10. 

The  ability  of  the  entire  procedure  derived  for  Schmidt's  statistic  in  estimating  the 
regression  coefficients  will  be  determined  for  a  best  and  worst  case  scenario.  The  best 
case  will  be  for  t2  =  .01  and  T  =  100,  while  the  worst  case  will  be  for  T2  =  10  and  T 
=  20.  Using  the  regression  model 

Yt  =  pxt  +  £t  (eqn  1.28) 

where 

ct  =  Pl£t-1  +  Pt-4  +  V  (efin  1-29> 
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for  Xt  =  .1  and  P  =  1  the  Yt  s  are  generated  for  various  values  of  pj  and  p_j.  Once 
the  Yt's  are  determined,  an  estimate  of  P  will  be  calculated.  Mean  square  error  will  be 
used  as  a  measure  to  compare  estimates  of  P  with  their  known  value  of  one.  This 
measure  will  indicate  the  efficiency  of  regression  parameters  estimated  by  the  procedure 
as  based  on  Schmidt's  statistic. 

E.  THESIS  ORGANIZATION 

The  thesis  is  organized  into  four  chapters.  Chapter  II  will  introduce  background 
theory'  utilized  to  derive  the  AR(1,4)  procedure.  In  particular,  the  Durbin- Watson  test 
for  the  AR(1)  process  will  be  addressed  along  with  the  Wallis  and  Schmidt  tests  for 
special  cases  of  the  AR(4)  process.  In  addition  to  test  procedures,  the  transformation 
for  autocorrelated  data  will  be  discussed. 

Chapter  III  will  initially  develop  the  test  for  the  AR(1,4)  process.  Specifically, 
upper  and  lower  bounds  will  be  developed  for  the  Boger  and  Schmidt  statistics.  After 
deriving  the  test,  the  data  transformation  for  the  AR(1,4)  process  will  be  addressed. 
Next,  results  will  be  presented  concerning  the  effectiveness  of  the  derived  results  for  the 
AR(1,4)  process.  Finally,  conclusions  and  areas  of  future  study  will  be  discussed. 
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II.  BACKGROUND  THEORY 


A.  GENERAL 

In  this  section  certain  concepts  underlying  all  time  series  processes  will  be 
defined.  These  definitions  will  lay  the  groundwork  which  will  lead  into  discussion  of 
specific  test  procedures.  In  particular,  the  Durbin-Watson,  Wallis  and  Schmidt  tests 
shall  be  discussed.  Discussion  of  these  procedures  will  act  as  a  prelude  to  the  next 
chapter's  development  of  a  test  and  correction  procedure  utilizing  an  extension  of 
Schmidt's  statistic. 

A  sample  of  error  residuals  ej,....ey,  where  the  index  denotes  a  point  in  time,  is 
called  a  time  series.  Each  observation  in  this  sample  is  a  realization  of  a  random 
variable  and  as  a  result  this  sequence  is  considered  to  be  a  discrete  stochastic  process. 
The  sequence  is  considered  discrete  since  observations  are  made  at  a  fixed  interval  of 
time.  Additionally,  ej,....ey  is  considered  as  a  finite  subsequence  of  an  infinite  series 
whose  joint  distribution  is  defined  by  these  finite  subsequences. 

Only  stationary  stochastic  processes  will  be  considered  in  this  paper.  A  process  is 
considered  stationary  if  the  joint  probability  distribution  of  the  residuals  ej,....ey  is 
unaffected  by  shifting  the  time  origin  forward  or  backward  by  k  time  units.  That  is, 
ej,....ey  and  e^+  i>—<ek  +  T’  are  identically  distributed  [Ref.  12). 

The  general  linear  model  with  autoregressive  disturbance  of  order  p  is  defined  as 

Yt  =  Xt'p  +  ct  (eqn  2.1) 

with 

ct  =  6lct-l  +  -  +  6pct-p  +  ^t  (ecl n  2-2) 

where  Yt  is  the  value  of  the  dependent  variable  at  time  t,  Xt'  is  a  (1  x  K)  vector 
representing  an  observation  on  K  non-stochastic  variables  at  time  t  and  P  is  a  (K  x  1) 
vector  of  coefficients  to  be  estimated.  Equation  2.2  implies  that  error  disturbances  in 
the  current  time  period  depend  on  those  in  previous  time  periods  and  another  error 
term  which  is  independent  with  zero  mean  and  constant  variance. 
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Autocovariance  defines  the  linear  relationship  between  members  of  a  stochastic 
sequence.  In  particular,  the  autocovariance,  y^,  between  ct  and  +  ^  is  given  by: 

Yk  =  E((£t'E^ct)Xct+k’E(£t  +  k)))  k  =  0,1,2,...  (eqn  2.3) 

where  E(£t)  is  the  mean  of  the  stochastic  process,  which  is  zero  in  our  application.  For 
a  stationary7  process  the  covariance  between  £t  and  £t+k  depends  on  k  and  not  on  the 
particular  point  in  time.  The  above  formula  completely  describes  the  autocovariance 
structure  of  the  stochastic  sequence  as  defined  above  [Ref.  13:  p.226]. 

Autocovariance  is  dependent  on  the  unit  of  measurement  for  the  underlying 
variable.  To  account  for  this,  the  y^'s  are  normalized  by  dividing  by  y*,  which  is  the 
variance  of  ev  Dividing  y^  by  */q  results  in  the  autocorrelation  function  of  £t: 

Pk  =  Yk/Yo>  k=  0,1,2,..  (eqn  2.4) 

From  the  above  definition  of  the  autocorrelation  function  it  is  obvious  that 

p0  =  1.  (eqn  2.5) 

In  order  to  gain  estimates  for  the  0's  or  autoregressive  parameters  in  terms  of  the 
autocorrelation  coefficients  or  p's,  the  Yule- Walker  equation  is  used  [Ref.  12:  p.  55). 
This  equation  for  the  above  general  linear  model  is  derived  by  left  multiplying  the 
equation 

£t  =  ei£t.i  +  ...  +  ep£t.p+nt  (eqn  2.6) 

by  £,_k>  and  then  taking  expectations  and  dividing  by  the  variance  Jq  of  £r  This 
results  in  the  Yule- Walker  equation 

Pk  =  eiPk-l  +  -  +  0pPk-p  fork*  1,2 .  (eqn  2.7) 

By  substituting  k=  1,2, ...p  into  the  above  equation,  a  set  of  linear  equations  for 
0j,02-  -0^  in  terms  of  Pj^-.-Pp  are  obtained.  By  substituting  estimates  for  P|.p->...Pp 
into  the  above  equations,  the  autoregressive  parameters  are  estimated.  Various  other 
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methods  exist  to  estimate  the  p's,  and  some  will  be  discussed  later.  This  result  is 
important  in  the  specification  of  the  autoregressive  process  and  will  help  in  the 
development  of  a  transformation  matrix  for  the  autocorrelated  data  [Ref.  12]. 

B.  DURBIN- WATSON  TEST 

This  section  summarizes  the  most  important  features  of  the  theory  for  the 
Durbin-Watson  test  (Ref.  2,14,15].  This  theory  is  the  basis  for  later  generalizations  of 
tests  for  autoregressive  disturbances. 

The  Durbin-Watson  test  is  based  on  the  statistic 


d  =  I(eret-i)2/Ict2=  e'Ae/e'e 

<*y 

where  e  =  y  -  XP  is  a  vector  of  least  squares  residuals  and  A  equals: 


(eqn  2.8) 
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For  certain  cases  where  regression  vectors  are  eigenvectors  of  the  matrix  A 
occurring  in  the  residuals  distribution,  the  statistic 


e'Ae  /  e’e  (eqn  2.9) 

provides  a  test  which  is  uniformly  most  powerful  against  certain  alternative  hypotheses 
[Ref.  16:  p.  88]. 

The  general  linear  model  was  previously  defined  as 
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Y  =  XP  +  £, 


(eqn  2.10) 


where  Y  is  a  (T  x  1)  matrix,  X  is  a  (T  x  K)  matrix  of  T  observations  on  K  variables 
and  c  is  a  (T  x  1)  matrix  of  errors.  The  least  squares  estimate  of  p  is  P  which  is  given 
by: 

(X'X^X'y.  (eqn  2.11) 

Also,  e,  the  vector  of  residuals  from  the  regression,  is  defined  as 

e  =  y  -  X$  =  My  (eqn  2.12) 

where  I  is  an  identity  matrix  of  order  T,  and 

M  =  I  -  X(X'X)‘1X'.  (eqn  2.13) 

It  can  be  verified  that  the  matrix  M  is  idempotent.  Substituting  XP  +  e  for  y  into  the 

formula  for  e  leads  to  the  result 

e  =  Me.  (eqn  2.14) 

This  result  implies  that 

d  =  e'MAMe  /  e'Me.  (eqn  2.15) 

The  distribution  of  the  above  statistic  is  dependent  on  the  (T  -  K)  non-zero  eigenvalues 
of  AM  which  are  denoted  by  S-.  Values  for  Sj  are  dependent  on  M  through  the  matrix 
X. 

The  matrix  A  is  a  (T  x  T)  symmetric  matrix  with  (T  -  1)  positive  eigenvalues 
which  for  T  greater  than  or  equal  to  three  lie  in  the  closed  interval  from  zero  to  four. 
The  eigenvalues  of  A  will  be  denoted  by  Xj  and  are  defined  as 

X-  =  2(1  -  cos(rcj/T))  j  =  1 . T-l.  (eqn  2.16) 
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Durbin  and  Watson  showed  that  if  the  et  are  assumed  to  be  normally  distributed, 
there  exists  an  orthogonal  linear  transformation  of  e  to  v  such  that 

d  =  Z6ivi2  /  Evi2<  (eqn  2.17) 

tz/  i»i 

where  Vj  are  independent  and  identically  distributed  normal  variables  with  zero  mean 
and  variance  <x2.  The  distribution  of  d  thus  depends  on  the  which  are  the 
eigenvalues  of  M  [Ref.  2:  p.  412]. 

Using  the  above  result  to  obtain  the  distribution  of  d  is  tedious,  since  the  8j's 
depend  on  the  X  matrix  and  the  statistic  would  have  to  be  calculated  for  each  new 
matrix.  Consequently,  if  it  is  desired  to  have  a  test  which  is  not  restricted  to  a 
particular  X  matrix,  the  distribution  would  have  to  be  determined  for  all  X  matrices. 
This  task  is  impossible  since  there  are  an  infinite  number  of  matrices. 

Durbin  and  Watson  avoided  this  problem  by  constructing  a  bounds  test  where 
the  upper  and  lower  bounds,  dU  and  d^,  respectively,  are  independent  of  the  particular 
X  matrix.  They  were  able  to  accomplish  this  by  determining  inequalities  on  the  in 
terms  of  the  eigenvalues  Xj  of  the  matrix  A,  where  the  8-'s  and  Xj's  have  been  arranged 
in  increasing  order.  This  result  implies  that 

dL<d<dy  (eqn  2.18) 

where 

dU  =  X  Xi+  Kv  i2  /  I*v  i2>  (ecln  2-19> 

(•/  Ut 

and 

T-K  .  r--K 

dL  =  IVi  /  X  V  i  •  ^ecln  2-20> 

/X  /  * 

Durbin  and  Watson  were  unable  to  find  the  exact  distribution  of  these  bounds  but 
approximated  the  distribution  with  the  first  four  terms  of  a  series  expansion  in  terms  of 
Jacobi  polynomials  using  a  Beta  distribution  as  a  weight  function. 

Using  the  statistic 

t  -  r  0 

X(  et-et-i  >  /  Set  ’  (ecin  2-21) 

-t-/ 
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and  asymptotic  results  it  can  be  shown  that 


d  =  2(  1  -  p  ), 


(eqn  2.22) 


where  p  is  an  estimate  of  the  autocorrelation  coefficient.  This  result  shows  that  when 
p  equals  zero,  indicating  no  autocorrelation,  d  equals  two.  In  addition,  when  p  equals 
one  indicating  positive  autocorrelation,  d  equals  zero  and  when  p  equals  negative  one, 
d  equals  four.  In  testing  the  null  hypothesis  of  no  autocorrelation  against  the 
alternative  hypothesis  of  positive  autocorrelation,  the  null  hypothesis  would  be  rejected 
if  the  calculated  value  for  d  is  less  than  the  lower  bound  calculated  and  accepted  if  d  is 
greater  than  the  upper  bound.  If  the  calculated  value  of  d  falls  between  the  upper  and 
lower  bounds  the  test  is  considered  inconclusive  and  there  is  insufficient  evidence  to 
accept  or  reject  the  null  hypothesis. 

Since  Durbin  and  Watson  first  calculated  upper  and  lower  bounds  for  their 
statistic  using  an  approximate  distribution,  a  method  has  been  developed  by  Imhoff 
which  allows  the  calculation  of  an  exact  distribution  for  a  quadratic  form  in  normal 
variables  [Ref.  7:  p.  419|.  To  briefly  indicate  how  the  Imhoff  distribution  is  utilized  in 
the  calculation  of  the  statistic's  bounds,  only  the  lower  bound  will  be  considered.  The 
calculation  of  the  upper  bound  can  be  done  in  an  analogous  manner.  In  order  to 
calculate  the  lower  bound,  the  eigenvalues  of  the  A  matrix,  X-,  are  calculated.  For  the 
A  matrix,  there  will  be  (T  -  1)  non-zero  eigenvalues  where  T  is  the  order  of  the  A 
matrix  and  the  one  zero  term  is  the  eigenvalue  which  corresponds  to  the  constant  term 
of  the  regression.  From  these  eigenvalues,  the  (T  -  K)  smallest  are  selected  to  compute 
the  statistic's  lower  bound  [Ref.  8:  p.70[.  The  value  for  K  is  equal  to  the  number  of 
regressors  including  a  constant  term.  The  basis  for  the  above  procedure  is  the  result 
proved  by  Durbin  and  Watson  which  determined  inequalities  of  the  eigenvalues  Sj  of 
the  matrix  AM  in  terms  of  the  eigenvalues  X-  of  the  A  matrix  or 


i  =  1,2,...T-K  , 


(eqn  2.23) 


where 


K'  =  K  -  1 


(eqn  2.24) 
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and  K'  is  the  number  of  regressors  not  including  a  constant  term.  Analogously,  to 
calculate  the  upper  bound  the  (T  -  K)  largest  eigenvalues  are  selected.  Adapting  the 


Imhoff  distribution  to  this  problem  leads  to 

F(dL)  =  1/2  -  1/rcj  ((  sin  <p(u))  /  (uto(u)))  du 

0 

(eqn  2.25) 

where 


7-K  . 

<p(u)  =  .5  J  tan'TUjU) 

/«/ 

(eqn  2.26) 

w(u)  =  ri(i  +  u:^u^)^ 

/!/ 

(eqn  2.27) 

and 


uj  =  Xj  -  dL  (i  =  1 . T-K) 

(eqn  2.28) 

Since  the  Xj's  are  known,  a  value  for  d^  is  selected  and  the  uj's  are  calculated. 
Using  the  Uj's,  the  integral  is  then  evaluated  numerically.  Numerical  integration  of  this 
infinite  integral  is  possible  since  Imhoff  has  proved  that  the  integral's  limit  as  u 
approaches  zero  is 


T-K 

l/2£  Uj. 

/*/ 

(eqn  2.29) 

Since  the  integral's  range  is  infinite,  there  will  be  a  truncation  error  associated  with  this 
integration.  Imhoff  has  shown  that  the  truncation  error  caused  by  integrating  over  the 
finite  range  from  0  to  u  inclusive,  can  be  held  to  any  arbitrary  level  p  by  taking 


u  =  (((T  -  K)/2)  rcpnVjl1/^  )g 
/-/ 

(eqn  2.30) 

where  g  =  2/(T-K).  Additionally,  the  numerical  integration  will  be  subject  to  error 
related  to  the  method  of  integration.  This  error  is  controlled  by  setting  the  error 
tolerance  within  the  program  used  to  integrate  the  function. 
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C.  EXTENSIONS  TO  HIGHER  ORDER  PROCESSES 

Durbin  and  Watson's  original  work  has  been  extended  to  higher  order 
autoregressive  processes.  Specifically,  Wallis  has  extended  it  to  a  special  case  of  an 
AR(4)  process  and  Schmidt  has  extended  it  to  a  second  order  process.  For  the 
regression  model 

Y  =  XP  +  £,  (eqn  2.31) 

with  an  error  specification 

ct  =  P£t-4  +  Ht»  (eqn  2-32) 

Wallis  generated  bounds  for  the  statistic 

r  r 

d4  =  X(et  *  et-4^  /Xet  =  e  A4e  /  e'e-  (ecln  2-33) 

Wallis'  derivation  of  this  result  is  analogous  to  Durbin  and  Watson's.  Using 
ImhofTs  distribution,  Wallis  was  able  to  calculate  bounds  for  this  special  case  of  the 
AR(4)  process  [Ref.  5:  p.  617],  Vinod,  in  a  manner  similar  to  Wallis,  further  extended 
these  results  to  include  bounds  for  the  same  special  cases  of  the  AR(2)  and  AR(3) 
processes  [Ref.  17]. 

In  order  to  account  for  an  AR(1)  or  AR(4)  disturbance,  the  data  must  be 
transformed.  In  this  regard,  for  the  AR(1)  process  the  dependent  variables  are 


transformed  as 

Yt*  =  YjO-pj2)1/2,  t=l  (eqn  2.34) 

and 

Y*  =  Yt-PlYM,  t=  2,3 . T.  (eqn  2.35) 

Similarly,  the  independent  variables  are  transformed  as: 

Xt*  =  Xt(l  -  Pj2)1/2,  t=l  (eqn  2.36) 
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and 


Xt*  =  Xt  -  pjXt_!,  t  =  2,3, . T.  (eqn  2.37) 

In  order  to  utilize  this  transformation,  an  estimate  for  pj  is  required.  Several 
methods  are  available  to  estimate  this  parameter.  The  method  of  choice  due  to  its 
favorable  properties  and  ease  of  computation  [Ref.  3]  is 

Pj  =  1  -  .5dj.  (eqn  2.38) 

For  the  special  case  of  the  AR(4)  process  the  transformation  is  similar. 
Specifically,  the  dependent  variables  are  transformed  as: 


Yt*  =  Yt(l  -  p/)1/2  t=  1,2,3, 4  (eqn  2.39) 

and 

Yt*  =  Yt  -  p4Yt_4  t  =  5 . T.  (eqn  2.40) 

The  independent  variables,  similarly  are  transformed  as: 

Xt*  =  Xt(l  -  p/)1/2  t=  1,2, 3, 4  (eqn  2.41) 

and 

Xt*  =  xt  -  p4Xt.4  t  =  5, T.  (eqn  2.42) 


As  in  the  AR(1)  case,  an  estimate  of  p4  is  required.  The  formula  for  this  estimate  is 
the  same  as  before  except  that  the  Wallis  statistic  is  used  in  lieu  of  the  Durbin-Watson. 
Once  either  of  these  transformations  is  utilized  after  the  detection  of  its  respective 
disturbance,  estimates  of  regression  parameters  can  be  expected  to  be  efficient  and 
unbiased. 

Schmidt  [Ref.  9]  generalized  Durbin  and  Watson's  preliminary  work  to  the  case 
of  a  second  order  autoregressive  process  whose  disturbance  is  of  the  form: 
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£t  =  Pl£t-1  +  f>2£t-2  +  Ht- 


(eqn  2.43) 


The  statistic  which  Schmidt  derived  to  detect  this  disturbance  was: 


d2  =  di  +  62 

where 


=  2(«t 

i*L 


-  e 


t- 


2 

t 


and 


s2  =  2(et  -  et-2>2  /  Set2' 

In  matrix  notation  this  statistic  is  represented  as 


(eqn  2.44) 


(eqn  2.45) 


(eqn  2.46) 


d2  =  e'(Aj  +  A2)e  /  e'e  (eqn  2.47) 

where  Aj  is  the  matrix  as  previously  defined  for  the  Durbin-Watson  statistic  and  A2  is 
defined  as 
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Following  results  proved  by  Durbin  and  Watson,  Schmidt  generated  the  statistic,  d2, 
utilizing  ImhofFs  distribution  and  the  eigenvalues  from  the  matrix  (Aj  +  A2).  In  a 
manner  similar  to  previous  developments,  he  further  showed  that 


<^2  “  2(2  -  Pj  -  P2)- 


(eqn  2.48) 


D.  SEQUENTIAL  TESTING 

In  development  of  a  generalization  of  the  Durbin-Watson  statistic,  Vinod, 
derived  a  scheme  of  sequential  tests  for  the  presence  of  autocorrelation  [Ref.  17].  In 
this  procedure  a  sequence  of  tests  is  used  to  determine  if  AR  processes  of  sequentially 
increasing  order  exist.  As  an  example,  the  null  hypotheses  for  the  four  tests  required 
to  detect  a  fourth  order  process  would  be: 

•  H0  :  pj  =  0, 

•  Hq  :  p2  =  0  given  p}  =  0, 

•  Hq  :  p3  =  0  given  Pj  =  P2  =  0, 

•  H0  :  p4  =  0  given  Pj  =  p2  =  P3  =  0. 

Thus  in  order  to  test  for  an  autoregressive  disturbance,  the  null  hypotheses  for 
sequentially  greater  orders  are  tested  until  one  of  the  null  hypotheses  is  rejected. 
Specifically,  at  the  step,  j  >  1,  the  hypothesis,  Hq  :  Pj  =  0  given  that  Pj  =  p2  = 
...  =  Pj_j  =  0  is  tested  against  the  two-sided  alternative  Ha  :  Pj  >  0  and  Pj  <  0. 
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III.  THE  AR(1,4)  PROCESS 


A.  INTRODUCTORY  ANALYSIS 

In  this  section  stationary  conditions  and  the  development  of  the  error  term  for 
the  AR(1,4)  process  will  be  presented.  If  stationary  conditions  do  not  hold  the 
AR(1,4)  model  is  invalid  and  a  model  which  considers  a  dynamic  specification  must  be 
pursued.  In  particular,  stationary  conditions  must  be  determined  for  the  regression 
model 


Yt=Xt'P  +  ct,  (eqn  3.1) 

whose  error  term  is  given  by: 

®  l^t- 1  ^4^t-4  ^t’  (eqn  3.2) 

As  discussed  in  previous  chapters,  the  Tjt's  are  assumed  to  be  distributed  normally  with 
zero  mean  and  constant  variance. 

Wise  [Ref.  18]  derived  a  method  to  determine  the  stationary  conditions  for  a 
stochastic  process  of  the  autoregressive  and  moving  average  types.  Following  his 
procedure,  it  is  possible  to  show  that  the  stationarity  conditions  for  a  fourth  order 
autoregressive  process  are 

•  2  +  pj  -  pj  +  2p4  >  0, 

•  3  +  p2  -  3p4  >  0, 

•  2  -  Oj  +  pj  +  2p4  >  0, 

•  1  -  Pi  -  P2  -  P3  -  P4  >  °- 

•  1  +  Pi  -  P2  +  P3  -  P4  >  °- 

Since  P2  and  p^  are  assumed  to  be  zero  in  the  AR(1,4)  process,  these  stationarity 
conditions  reduce  to 

•  2  +  Pj  +  2p4  >  0, 

•  3(1  -  p4)  >  0, 

•  2  -  pj  +  2p4  >  0, 

•  1  -  Pi  -  P4  >  0, 

•  1  +  Pi  -  Pq  >  0. 
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The  error  process  for  a  general  fourth  order  process  is  defined  by  the  relation: 


Et  =  9lEt-l+92Et-2  +  93Et-3  +  94Et-4-  (e<In  3l3) 

If  the  relationship  is  multiplied  through  by  and  expectations  are  taken,  the 


difference  equation 

?k =  9 1  Yk- 1  +  02?k-2  +  -  +  64Yk-4  k  >  0 

(eqn  3.4) 

results,  where  y  was  previously  defined  as  the  autocovariance.  If  this  equation  is 


divided  through  by  Jq,  the  result  is 

Pk  =  9 1  Pk- 1  +  02Pk-2  +  • +  04Pk-4- 

(eqn  3.5) 

which  is  the  general  form  of  the  Yule-Walker  equation  and  defines  the  autoregressive 
parameters,  0's,  in  terms  of  the  autocorrelations,  p's.  If  k  =  1,  2,  3,  4,  is  substituted 
into  this  equation  a  set  of  linear  equations  is  obtained  for  0j,  ©2,  03  and  04  in  terms  of 
Pj,  p->,  p3  and  p4.  Specifically,  these  equations  are 


Pj  =  0j  +02pj  +  03p2  +  04p3, 

(eqn  3.6) 

P2=0lPl+02  +  e3Pl+e4P2- 

(eqn  3.7) 

P3  =  9lP2  +  92Pl  +93  +  94Pl> 

(eqn  3.8) 

and 


P4=9lP3  +  92P2  +  93Pl  +  94- 

(eqn  3.9) 

Since  the  error  process  for  the  AR(1,4)  process  is  specified  as 

Et  =  9lEt-l  +04ct-4+nt>  (eqn  3.10) 
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where  ©2  and  63  are  assumed  equal  to  zero,  the  above  equations  can  be  reduced  to 


Pj  =  0j  +  04P3,  (eqn3.11) 

P2  =  6 1 P 1  +  6qP2’  (eqn  3.12) 

P3  =  ©  1P2  +  ®qPl>  (eqn  3.13) 

and 

P4  =  0 1P3  +  64-  (eqn  3.14) 

Further  since  it  is  assumed  that  P[  and  P3  are  not  present  in  the  disturbance,  the  above 
equations  can  be  finally  reduced  to 

Pj  =  0j  (eqn  3.15) 

and 

P4  =  64-  (eqn  3.16) 

These  final  equations  allow  the  AR(1,4)  error  term  to  be  represented  as: 

"b  Pq£j_4  "l"  (eqn  3,17) 

B.  BOGER'S  STATISTIC 

Boger  proposed  a  statistic  which  was  a  variation  of  the  Durbin-Watson  statistic 
to  test  for  the  AR(1,4)  disturbance.  This  statistic  is  of  the  form 

d  1 4  =  E  (et-et.  j  -et.4)2  /  £et2.  (eqn  3. 1 8) 

*•5  t1*/ 
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In  matrix  notation  this  statistic  is  denoted  as 


dl,4  =  e  Al,4e  /  e'e 
where  A  j  4  equals 


(eqn  3.19) 
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Bounds  for  Boger's  statistic  could  not  be  determined  because  its  A  matrix  as 
defined  above  did  not  satisfy  properties  as  established  by  Durbin  and  Watson.  In 
particular,  Durbin  and  Watson  [Ref.  2]  showed  that  the  matrix  product  MA,  where  A 
is  any  real  symmetric  matrix  has  (T-K)  real  positive  eigenvalues  of  which  K  are  zero 
eigenvalues.  T  and  K  are  the  dimensions  of  the  matrix  X.  Specifically,  Durbin  and 
Watson  showed  for  their  statistic  whose  A  matrix  equals 
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that  there  were  T-l  positive  eigenvalues.  As  previously  discussed,  they  also  showed 
that  if  the  positive  eigenvalues  from  both  MA  and  A  are  placed  in  increasing  order, 
bounds  for  statistics  of  the  form 


d  =  e'Ae  /  e'e, 


(eqn  3.20) 


can  be  calculated. 

Where  Boger's  statistic  violated  these  properties  was  that  for  all  values  of  K  or 
the  number  of  independent  variables,  there  were  not  sufficient  eigenvalues  in  the  A 
matrix  to  bound  those  eigenvalues  in  the  matrix  MA.  Specifically,  the  A  matrix  for 
Boger's  statistic  had  one  zero  and  three  negative  eigenvalues  where  Durbin  and 
Watson's  had  one  zero  eigenvalue.  Since  the  matrix  MA  had  T-K  positive  eigenvalues 
and  Boger's  A  matrix  only  had  T-4  positive  eigenvalues  instead  of  T-l,  when  the 
number  of  independent  variables  was  less  than  or  equal  to  three,  there  were  not 
sufficient  eigenvalues  from  Boger's  A  matrix  to  bound  those  from  the  matrix  MA. 
Since  this  violated  Durbin  and  Watson's  theory,  upper  and  lower  bounds  could  not  be 
calculated  for  Boger's  statistic,  when  the  number  of  independent  variables  was  less  than 
three. 

C.  SCHMIDT'S  STATISTIC  EXTENDED  TO  THE  AR(1,4)  PROCESS 

Since  it  wras  desirable  to  have  a  statistic  which  was  not  limited  by  the  number  of 
independent  variables,  a  variation  of  Schmidt's  statistic  to  fit  the  AR(1,4)  process  w'as 
investigated.  Schmidt  was  able  to  show  that  the  statistic  d^  where 

d2  =  dj  +  <>2  (eqn  3.21) 


and 


(eqn  3.22) 


provides  a  test  to  detect  a  second  order  autoregressive  disturbance. 

Using  Schmidt's  research  this  paper  proposes  that  it  is  possible  to  detect  the 
AR(1,4)  process  by  way  of  the  statistic  dj  4  where: 


d  l  ,4  =  d!  +  d4 


(eqn  3.23) 
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or 


In  matrix  notation  this  statistic  is  given  as 


(eqn  3.24) 


e'A14e  /  e'e 


(eqn  3.25) 


where 


Al,4  -  Aj  +  A4 


(eqn  3.26) 


or  A  j  4  equals: 
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Results  established  by  ImhofF  and  Durbin  and  Watson  were  used  to  determine 
upper  and  lower  bounds  for  the  statistic  dj  4.  Koerts  and  Abrahamse  [Ref.  8]  have 
provided  Fortran  66  code  which  allows  the  determination  of  the  distribution  of  a 
quadratic  form  in  normal  variables  by  way  of  ImhofFs  distribution.  This  code  was 
updated  to  Fortran  77  standards  and  modified  to  allow  the  calculation  of  the  dj  4 
statistic.  The  algorithm  utilizes  ImhofFs  distribution  to  determine  five  percent 
significance  points  of  the  statistic  by  successively  halving  the  range  of  the  function 
until 


F(dL)  =  .05, 


(eqn  3.27) 
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where  F(dj^)  is  the  Imhoff  distribution  function  as  defined  in  the  previous  chapter. 
Truncation  and  integration  errors  are  controlled  by  input  parameters.  It  was  found 
that  if  these  errors  were  set  to  less  than  .0001,  the  program  required  inordinate 
amounts  of  CPU  time.  At  a  .0001  level  of  accuracy,  the  five  percent  significance  points 
utilized  S2500  of  computer  resources.  At  a  .00001  level  of  accuracy,  approximately  a 
six  fold  increase  in  resources  can  be  expected.  Five  percent  significance  points  for  the 
statistic  dj  ^  are  listed  in  Table  I.  In  Table  I,  K  is  the  number  of  regressors  including 
the  intercept  term  and  T  is  the  number  of  observations. 

Currently,  if  an  AR(1,4)  process  is  suspected  to  exist  in  the  data  a  two-step 
procedure  is  generally  used.  Specifically,  this  involves  first  testing  and  correcting  for 
the  AR(1)  component  of  the  disturbance  and  then  doing  the  same  for  the  AR(4) 
component.  The  AR(1,4)  procedure  proposes  to  test  and  correct  the  AR(1,4) 
disturbance  in  one  step  where  both  components  are  simultaneously  corrected.  Once 
autocorrelation  has  been  determined  to  exist  within  the  data  for  both  procedures, 
values  for  the  autocorrelation  coefficients  Pj  and  p4  must  be  determined.  In  the 
two-step  procedure,  Pj  is  calculated  by  way  of  the  formula 


Pj  —  1  -  .5dj, 


(eqn  3.28) 


where  dj  is  the  Durbin- Watson  statistic.  After  an  appropriate  transformation,  the  data 
are  tested  in  a  sequential  manner  as  discussed  in  the  previous  chapter  until  a  higher 
order  disturbance  is  detected.  If  a  fourth  order  disturbance  is  detected,  p4  is  estimated 
by  the  formula 


(eqn  3.29) 


p4  =  1  -  .5d4 


where  d4  is  the  statistic  tabulated  by  Wallis.  As  can  be  seen  by  this  procedure,  the 
calculation  of  values  for  p  j  and  p4  are  in  a  sense  conditional  upon  each  other.  That  is 
since  pj's  value  is  used  as  part  of  the  original  data  transformation  and  p4  is  determined 
by  way  of  the  Wallis  statistic  which  is  calculated  from  the  transformed  data,  the  value 
of  p4  depends  on  the  value  of  pj. 

There  are  two  important  assumptions  concerning  the  two-step  procedure  as  it 
relates  to  the  AR(l,4)  disturbance.  The  first  of  these  is  that  the  Durbin-Watson 
statistic  will  be  used  to  detect  the  AR(1,4)  process.  That  is,  that  an  individual  who  is 
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TABLE  I 

TABLE  OF  BOUNDS  FOR  THE  SCHMIDT  STATISTIC 


K  -- 


4 


T 

L 

U 

L 

U 

L 

U 

15 

2.38 

2.84 

2.19 

3.05 

2.03 

3.29 

16 

2.43 

2.87 

2.26 

3.06 

2.10 

3.30 

17 

2.48 

2.91 

2.31 

3.10 

2.17 

3.30 

18 

2.53 

2.94 

2.36 

3.13 

2.23 

3.30 

19 

2.57 

2.97 

2.42 

3.15 

2.27 

3.32 

20 

2.61 

3.00 

2.47 

3.16 

2.33 

3.34 

21 

2.65 

3.02 

2.51 

3.18 

2.37 

3.35 

22 

2.68 

3.03 

2.55 

3.20 

2.42 

3.35 

23 

2.72 

3.05 

2.58 

3.22 

2.46 

3.36 

24 

2.75 

3.08 

2.62 

3.24 

2.50 

3.37 

25 

2.78 

3.09 

2.65 

3.26 

2.54 

3.38 

26 

2.80 

3.11 

2.69 

3.26 

2.57 

3.39 

27 

2.82 

3.12 

2.71 

3.28 

2.60 

3.40 

28 

2.85 

3.14 

2.74 

3.29 

2.63 

3.41 

29 

2.87 

3.15 

2.77 

3.30 

2.66 

3.42 

30 

2.90 

3.17 

2.80 

3.31 

2.69 

3.43 

31 

2.92 

3.17 

2.82' 

3.32 

2.72 

3.44 

32 

2.94 

3.19 

2.84 

3.33 

2.74 

3.44 

33 

2.96 

3.20 

2.86 

3.33 

2.77 

3.45 

34 

2.98 

3.21 

2.88 

3.34 

2.79 

3.46 

35 

2.99 

3.22 

2.90 

3.35 

2.81 

3.46 

36 

3.01 

3.23 

2.92 

3.36 

2.83 

3.47 

37 

3.03 

3.24 

2.94 

3.36 

2.85 

3.48 

38 

3.04 

3.25 

2.96 

3.37 

2.87 

3.48 

39 

3.05 

3.26 

2.97 

3.38 

2.89 

3.49 

40 

3.07 

3.27 

2.99 

3.38 

2.91 

3.49 

45 

3.13 

3.31 

3.05 

3.41 

2.98 

3.50 

50 

3.18 

3.34 

3.12 

3.43 

3.05 

3.52 

55 

3.23 

3.37 

3.16 

3.46 

3.10 

3.54 

60 

3.26 

3.40 

3.21 

3.47 

3.15 

3.55 

65 

3.30 

3.42 

3.24 

3.49 

3.19 

3.56 

70 

3.33 

3.44 

3.28 

3.51 

3.23 

3.57 

75 

3.35 

3.46 

3.31 

3.52 

3.26 

3.58 

80 

3.38 

3.47 

3.34 

3.53 

3.30 

3.59 

85 

3.40 

3.49 

3.36 

3.54 

3.32 

3.59 

90 

3.42 

3.50 

3.38 

3.55 

3.34 

3.60 

95 

3.43 

3.51 

3.40 

3.56 

3.36 

3.61 

100 

3.45 

3.53 

3.42 

3.57 

3.38 

3.62 

34 


TABLE  I 

TABLE  OF  BOUNDS  FOR  THE  SCHMIDT  STATISTIC  (CONT'D.) 


K--  5  6 


r 

L 

U 

L 

U 

15 

1.87 

3.47 

1.73 

3.74 

16 

1.95 

3.47 

1.82 

3.73 
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3.47 
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40 
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45 
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3.18 

3.63 
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3.69 

75 

3.22 

3.64 

3.18 

3.69 

80 
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3.64 

3.21 

3.70 

85 

3.28 

3.65 

3.24 

3.70 

90 

3.30 

3.66 

3.27 
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95 

3.33 

3.66 

3.29 

3.71 

100 

3.35 

3.66 

3.32 

3.71 

35 


testing  for  the  AR(1,4)  process  will  not  calculate  the  Wallis  statistic  if  the 
Durbin-Watson  statistic  fails  to  detect.  This  assumption  should  be  fairly  accurate  in 
practice  as  an  individual  would  most  likely  decide  that  autocorrelation  is  not  present  if 
the  Durbin-Watson  statistic  indicated  such.  The  second  assumption  which  is  closely 
related  to  the  first,  concerns  the  order  of  calculation  for  pj  and  p4.  Specifically,  as 
previously  shown,  the  formula  derived  from  the  Schmidt  statistic  contains  a  Pj  and  p4 
term.  Since  this  is  the  case,  different  results  might  have  been  obtained  if  p4  had  been 
calculated  by  way  of  the  Wallis  statistic  formula  first,  after  which  Pj  would  have  been 
calculated  by  the  Schmidt  statistic  formula.  This  possibility  w'as  not  investigated  in 
this  paper,  but  Pj  was  calculated  by  way  of  Durbin  and  Watson's  statistic  after  which 
p4  was  determined  from  the  Schmidt  formula. 

For  the  AR(1,4)  procedure  values  for  Pj  and  p4  wall  be  calculated  in  a  similar 
manner.  Specifically,  it  is  proposed  that  Pj  will  be  calculated  by  the  formula 

Pj  =  1  -  .5dj  (eqn  3.30) 

and  p4  by  the  formula 

p4  =  2  -  pj  -  .5dj  4.  (eqn  3.31) 

The  formula  for  p4  is  derived  from  the  original  definition  of  dj  4  or 


(eqn  3.32) 


If  the  formula  for  dj  4  is  factored  and  asymptotic  arguments  applied,  after  appropriate 
rearrangement  of  terms,  the  formula  for  p4  is  achieved.  The  validity  of  these  formulas 
is  argued  in  a  manner  similar  to  that  for  the  two-step  procedure.  The  difference  is  that 
values  for  pj  and  p4  are  determined  in  a  single  step.  An  advantage  to  this  method  is 
that  a  value  for  p4  is  calculated  prior  to  the  data  transformation,  thus  its  value  will  not 
be  influenced  by  any  peculiarities  which  may  have  occurred  during  the  data 
transformation. 

Once  estimates  for  pj  and  p4  have  been  determined,  a  method  to  utilize  them  in 
correcting  the  AR(1,4)  disturbance  must  nowr  be  determined.  As  previously  discussed, 
this  requires  that  a  matrix  P  be  found  such  that  P'P  =  ly'*  where  is  the  inverse  of 


36 


the  sample  covariance  matrix.  Once  the  P  matrix  is  determined,  the  data  can  be 
transformed  such  that 


Y*  =  PY,  (eqn  3.33) 

X*  =  PX,  (eqn  3.34) 

e  =  Pe,  (eqn  3.35) 

and  least  squares  applied  to  the  model 

Y  =  X  P  +  e  ,  (eqn  3.36) 

where  the  estimator  for  p  is 

£  =  (X'P'PX)'IX'P'PY.  (eqn  3.37) 


If  P  is  determined  from  the  transformed  data,  it  should  be  more  efficient  than  estimates 
from  ordinary  least  squares. 

Beach  and  MacKinnon  [Ref.  11],  in  their  paper  on  maximum  likelihood 
estimation  of  a  second  order  process,  defined  the  P  matrix  as 
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They  proposed  a  secondary  method  to  determine  values  for  a,  b  and  c  where  P'P  is 
constrained  to  be  proportional  to  the  inverse  of  the  variance-covariance  matrix  as 
given  be  Siddiqui  [Ref.  10]  and  a,  b  and  c  are  solved  for  in  the  implied  restrictions.  In 
a  manner  similar  to  that  proposed  by  Beach  and  MacKinnon,  this  paper  will  solve  for 
the  P  matrix  of  the  fourth  order  process. 

From  Siddiqui's  paper  on  the  inversion  of  the  sample  covariance  matrix  it  was 
possible  to  determine  that  the  inverse  of  the  sample  covariance  matrix,  V|/'*,  equals 
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where  a j  =  p j  and  a^  =  p^.  This  matrix  was  obtained  from  the  matrix  for  a  fourth 
order  process  after  setting  the  terms  for  a2  and  a^  equal  to  zero,  where  a^p-,  and 

a3  =  f*3* 

The  P  matrix  for  the  AR(1,4)  process  is  defined  such  that  P  equals: 
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a 


.  0 

be . 0 

def . 0 

ghij . 0 

-P4  0  0  -Pi  1 . 0 

iP . *P4  0  0  *Pj  1 

By  following  Siddiqui's  method  in  the  context  as  proposed  by  Beach  and  MacKinnon, 
10  equations  in  10  unknown  are  defined 

a2  +  b2  +  d2  -l-  g2  =  1  -  p42,  (eqn  3.38) 

be  +  de  +  gh  =  -p  i ,  (eqn  3.39) 

fd  +  ig  =  0,  (eqn  3.40) 

gj  =  -Pjpq,  (eqn  3.41) 

c2  +  e2  +  h2  =  1  +  pj2  -  p42,  (eqn  3.42) 

fe  +  hi  =  -p j,  (eqn  3.43) 
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jh  =  0, 


(eqn  3.44) 


f2  +  i2  =  1  +  p,2  -  p42, 

(eqn  3.45) 

ij  =  -Pi, 

(eqn  3.46) 

i2  =  1  *  P42 

(eqn  3.47) 

Solving  these  equations  for  the  10  unknowns,  the  elements  of  the  previously  defined  P 
matrix  are  determined  where 


a  =  V  (1  -  P42  -  b2  -  d2  -  g2), 

(eqn  3.48) 

b  =  (-pj  -  gh  -de)/  c, 

(eqn  3.49) 

c  =  V  (1  +  pj2  -  p42  -  h2  -  e2), 

(eqn  3.50) 

d  =  (-Pi2P4V(1  -  P42))/((l  -  P42)  V  ((1  -  P42)2  -  Pj2P42)), 

(eqn  3.51) 

e  =  (-PjV  (i  -  P42))  /  (V  ((1  -  P42)2  -  Pi2P42)), 

(eqn  3.52) 

f  =  V  (((1  -  P42)2  -  Pj2P42)  /  (1  -  P42)), 

(eqn  3.53) 
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g  =  -P1P4  /  J  (1  -  P42), 


(eqn  3.54) 


h  =  0, 


(eqn  3.55) 


i  =  -Pi  /  V  (1  -  Pq2), 


(eqn  3.56) 


i  =  V  (i  -P42). 


(eqn  3.57) 


The  previously  defined  matrix  with  substituted  values  as  calculated  above  satisfies  the 
required  property  that  P'P  = 

D.  EFFECTIVENESS  OF  THE  AR(1,4)  PROCEDURE 

As  with  any  new  procedure,  it  is  important  that  its  effectiveness  relative  to  older 
procedures  be  determined.  In  this  case  the  AR(1,4)  procedure  will  be  compared  to  the 
two-step  procedure.  Specifically,  these  procedures  shall  be  compared  in  three  areas: 

•  Accuracy  of  estimate  for  pq, 

•  Ability  to  detect  the  AR(1,4)  disturbance 

•  Accuracy  of  estimates  for  regression  parameters 

Two  simulation  programs  written  in  APL  were  utilized  to  generate  data  for  the 
analysis  required  to  make  the  comparison.  These  programs  are  displayed  in  the 
appendix.  Program  Stats  is  used  to  generate  data  for  the  comparison  of  the  estimates 
of  pq  and  the  ability  to  detect  the  AR(1,4)  disturbance.  This  program  utilizes  the  APL 
programs  Norrand  and  Unirand  for  the  generation  of  normally  distributed  random 
numbers.  Norrand  and  Unirand  are  generic  APL  random  number  generators  on  the 
IBM  370/3033 AP  computer  system  which  was  used  to  perform  the  simulations.  After 
the  program  generates  regression  residuals,  the  Durbin-Watson,  Wallis  and  Schmidt 
statistics  are  calculated.  The  Durbin-Watson  and  Schmidt  statistics  are  then  compared 
to  their  respective  bounds  for  the  selected  number  of  observations  and  independent 
variables.  If  the  statistic  detects  the  disturbance,  a  counter  will  be  incremented.  In 
addition  to  determining  whether  a  detection  is  made  or  not,  the  statistics  are  used  in 
previously  discussed  formulas  to  calculate  values  for  pj  and  pq. 
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Program  Regress  is  the  second  simulation  used.  Its  specific  purpose  is  to 
generate  data  which  will  allow  the  comparison  of  regression  parameters  from  both  the 
two-step  and  AR(1,4)  procedures.  This  program  as  does  the  previous,  utilizes  the 
Norrand  and  Unirand  routines  to  generate  normally  distributed  random  numbers.  In 
the  program,  values  for  one  independent  variable  and  intercept  term  were  generated  for 
a  predetermined  number  of  observations,  from  a  normal  distribution  with  zero  mean 
and  a  variance  of  .0625.  After  generating  the  error  term  from  the  AR(1,4)  model  for 
specified  values  of  Pj  and  p^,  the  dependent  variables  are  determined  by  way  of  the 
general  linear  model  with  a  value  for  P  of  one.  Using  this  data,  estimates  of  P  are 
determined  by  way  of  ordinary  least  squares,  the  two-step  procedure  and  the  AR(1,4) 
procedure. 

As  previously  discussed  in  Chapter  I,  it  was  decided  to  determine  how  well  the 
AR(1,4)  procedure  performs  relative  to  the  two-step  procedure  for  a  best  case  and 
worst  case  scenario.  For  the  best  case  scenario,  the  number  of  observations  was  100 
and  the  error  term  was  generated  from  a  normal  distribution  with  zero  mean  and  a 
variance  of  .01.  For  the  worst  case  the  number  of  observations  was  20  and  the 
variance  of  the  normal  distribution  was  10.  For  both  the  best  and  worst  case 
scenarios,  200  replications  were  performed  for  17  combinations  of  Pj  and  p^.  The 
number  of  combinations  which  could  be  investigated  was  limited  by  the  initial  terms  of 
the  P  matrix  which  contained  square  roots.  The  problem  was  that  for  certain 
combinations  of  pj  and  p^  the  term  under  the  square  root  became  negative.  The  17 
combinations  of  Pj  and  p^  which  were  investigated  were 
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14  P|  —  .7  p4  —  .3 

15  Pj  =  .7  p4  =  *5 

16  pj  =  .9  p4  =  .1 

17  pj  =  .9  p4  =  .3  . 

For  a  given  value  of  p  j ,  any  value  of  p4  greater  than  the  displayed  values  for  p4 
caused  the  problem  with  terms  under  square  roots  in  the  P  matrix.  In  subsequent 
graphs  used  to  display  results,  each  combination  of  Pj  and  p4  will  be  referred  to  by  the 
number  as  assigned  above.  Thus  combination  one  is  given  where  Pj  =  .1  and  p4  = 
.1.  The  effectiveness  of  the  procedure  in  handling  negative  autocorrelation  was  not 
investigated. 

Mean  square  error  (MSE)  was  selected  to  be  the  measure  by  which  the  best 
procedure  would  be  selected.  In  the  comparison  of^s  from  both  procedures,  MSE 
was  estimated  separately  for  the  slope  and  intercept  terms  utilizing  the  formula 


(eqn  3.58) 


where  u  was  the  actual  value  of  the  slope  or  intercept  and  n  was  the  number  of 
simulation  replications.  Specifically,  the  value  of  u  was  unity  for  this  application.  In  a 
similar  manner,  the  \fsfe  for  the  vector  p  consisting  of  the  slope  and  intercept  term 
was  determined  via  the  formula 

MSE  =  £  ((f0  -  l)(fj  -  l))2/n  (eqn  3.59) 

where  ^  is  the  intercept  term  and  Pj  is  the  slope  term.  Mean  square  error  for  the 
comparison  of  p4's  was  calculated  in  a  similar  manner. 

Power  curves  were  constructed  to  allow  the  determination  of  which  procedure 
could  best  detect  the  AR(1,4)  disturbance.  These  curves  were  constructed  by  plotting 
for  each  Pjand  p4  combination  the  number  of  times  each  procedure  detected  the 
disturbance  out  of  100  possible  attempts.  As  previously  mentioned,  if  a  trial  turned 
out  to  bt  inconclusive,  it  was  not  counted  as  a  detection.  The  number  of  times  for 
each  Pj.p4  combination  that  the  procedures  were  inconclusive  was  plotted  to  enable 
readers  who  consider  inconclusive  trials  as  detections  to  reevaluate  the  power  curves. 
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E.  RESULTS 

Comparisons  of  mean  square  error  for  estimates  of  are  shown  in  Figures  3.1 
and  3.2.  Figure  3.1  displays  mean  square  errors  for  the  worst  case  scenario  while 
Figure  3.2  displays  MSEs  for  the  best  case  situation.  Mean  squared  errors  for  both  the 
best  and  worst  case  scenario  are  consistent  for  both  the  AR(1,4)  and  two-step 
procedure  estimates.  For  the  best  case  scenario  the  two-step  estimate  is  clearly 
superior  for  Pj,p^  combinations  greater  than  nine.  MSEs  for  two-step  and  AR(1,4) 
procedure  estimators  for  the  worst  case  scenario  diverge  for  p  j  ,p^  combinations  greater 
than  15.  This  indicates  that  for  the  worst  case  scenario,  the  AR(1,4)  procedure 
estimate  for  p^  can  be  expected  to  perform  as  well  as  the  two-step  estimate. 

As  previously  discussed,  an  individual  using  the  two-step  procedure  to  test  for  an 
AR(1,4)  disturbance  might  initially  use  the  Durbin-Watson  statistic.  If  the  test 
indicated  that  autocorrelation  did  not  exist,  further  testing  with  the  Wallis  statistic 
might  not  be  considered.  Since  this  was  the  case  assumed,  power  curves  were  created 
which  compared  the  Schmidt  and  Durbin-Watson  statistics.  The  Schmidt  statistic  is 
considered  as  part  of  the  AR(1,4)  procedure  while  the  Durbin-Watson  is  part  of  the 
two-step  procedure.  For  the  best  case  scenario  with  T  =  100  and  x  =  .01,  shown  in 
Figure  3.4,  the  Schmidt  statistic  except  for  the  sixth  combination  is  more  powerful  than 
the  Durbin-Watson  statistic.  For  this  scenario,  both  statistics  are  equally  powerful  for 
Pj.Pq  combinations  greater  than  10.  For  the  worst  case  situation,  shown  in  Figure  3.3, 
both  statistics  are  equally  capable  across  all  combinations  in  detecting  the  AR(1,4) 
disturbance.  As  previously  discussed  the  power  curves  for  both  procedures  do  not 
include  inconclusive  test  results.  Figures  3.5  and  3.6  display  the  number  of  times  that 
the  test  statistics  were  inconclusive  for  the  worst  and  best  case  scenarios  respectively. 

Figures  3.7  through  3.12  show  comparisons  of  mean  square  errors  for  the 

regression  parameters  of  the  two-step  and  AR(1,4)  procedures.  Figures  3.7  and  3.8 

show  \ISEs  calculated  using  both  Pq  and  Pj  for  worst  and  best  case  situations, 

respectively.  For  the  worst  case  with  T  =  20,  x  =  10  scenario,  shown  in  Figure  3.7, 
A 

MSEs  for  both  procedures  follow  similar  trends  which  are  consistent  with  each  other. 
Figures  3.9  and  3.10  which  display  MSEs  for  the  intercept  and  slope  coefficients, 
respectively,  show  trends  which  support  the  previous  findings  of  Figure  3.7.  Mean 
square  errors  for  the  best  case  scenario  are  shown  in  Figures  3.8,  3.11  and  3.12.  As 
shown  in  Figure  3.8,  both  procedures  follow  similar  trends  for  Pj.P/j  combinations  up 
to  15.  Clearly  the  two-step  procedure  diverges  at  the  Pj.Pq  combination  15.  Figures 
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Figure  3.1  MSB  Comparison  of  Estimators  for  p_j  Worst  Case 


Figure  3.2  MSE  Comparison  of  Estimators  for  Best  Case 
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Figure  3.3  Power  of  Statistics  Worst  Case 


Figure  3.4  Power  of  Statistics  Best  Case 
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Figure  3.5  Inconclusiveness  of  Statistics  Worst  Case 


Figure  3.6  Inconclusiveness  of  Statistics  Best  Case 
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3.11  and  3.12  indicate  that  this  divergence  is  related  to  a  divergence  of  the  intercept 
estimate  for  the  two-step  procedure.  Overall,  it  appears  that  regression  parameters 
from  both  procedures  are  consistent  with  each  other  and  that  one  procedure  is  not 
superior  to  the  other  with  regards  to  the  estimation  of  regression  parameters. 

F.  CONCLUSIONS  AND  RECOMMENDATIONS 

Overall,  the  AR(1,4)  procedure  as  previously  developed  is  just  as  capable  as  the 
two-step  procedure  in  correcting  the  AR(1,4)  disturbance  for  the  Pj.p^j  combinations 
for  which  it  can  be  used.  The  Schmidt  statistic  as  part  of  this  procedure  is  effective  in 
detecting  the  AR(1,4)  disturbance.  In  certain  situations,  its  abilities  surpass  those  of 
the  Durbin-Watson  statistic  from  the  two-step  procedure.  In  particular,  its  capabilities 
arc  noteworthy  at  higher  values  of  T,  lower  values  of  and  lower  p  j  ,p^  combinations. 
It  appears  from  the  results  that  estimates  of  Pj  calculated  via  the  two-step  procedure 
are  better  than  estimates  from  the  AR(1,4)  procedure,  at  least  for  the  order  of 
calculation  of  p_j  which  was  tested.  If  an  individual  suspects  that  an  AR(1,4) 
disturbance  is  present  in  the  data  then  it  is  recommended  that  the  Schmidt, 
Durbin-Watson,  and  Wallis  statistics  be  calculated.  The  Schmidt  statistic  should  be 
utilized  to  detect  the  disturbance  while  the  Durbin-Watson  and  Wallis  statistics  are 
used  to  calculate  values  of  pj  and  p^.  If  the  values  for  Pj  and  p_j  are  acceptable,  then 
they  should  be  used  in  the  AR(1,4)  P  matrix  to  correct  for  the  disturbance.  If  the  pj 
and  p^j  values  fall  out  of  the  acceptable  range  for  the  AR(1,4)  procedure  then  the 
two-step  procedure  must  be  used. 

G.  AREAS  OF  FUTURE  RESEARCH 

Three  areas  of  future  research  are  suggested  as  a  result  of  this  research.  The  first 
of  these  would  be  to  find  a  P  matrix  for  the  AR(1,4)  procedure  which  was  not  limited 
by  values  of  pj  and  p^.  A  second  area  of  research  would  be  to  derive  a  correction 
procedure  for  a  general  fourth  order  autoregressive  process.  This  derivation  would 
follow  closely  what  was  done  above  for  the  AR(1,4)  process. 

A  final  area  of  research  would  be  to  investigate  the  maximum  likelihood 
estimation  of  a  general  AR(4)  or  an  AR(1,4)  process.  Beach  and  MacKinnon  [Ref.  11] 
have  investigated  the  second  order  process  and  suggest  possible  ways  to  extend  this 
procedure  to  the  general  fourth  order  case. 
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Figure  3.7  Comparison  of  Two-step  and  AR(1,4)  Procedure's  MSEs  Worst  Case 


✓N 

Figure  3.8  Comparison  of  Two-step  and  AR(1,4)  Procedure's  MSEs  Best  Case 


49 


Figure  3.9  Comparison  of  Two-step  and  AR(1,4)  Procedure's  Intercept  MSEs  Worst  Case 


A 

Figure  3.10  Comparison  of  Two-step  and  AR(1,4)  Procedure's  Slope  MSEs  Worst  Case 
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Figure  3.11  Comparison  of  Two-step  and  AR(1,4)  Procedure's  Intercept  MSEs  Best  Case 
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Figure  3.12  Comparison  of  Two-step  and  AR(1,4)  Procedure's  Slope  MSEs  Best  Case 
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APPENDIX 
PROGRAM  LISTINGS 


C 

c 


THIS  PROGRAM  COMPUTES  THE  LOWER  BOUND  FOR  THE  SCHMIDT  STATISTIC 
USING  IMHOFF'S  INTEGRAL.  THE  FILE  MAT  IS  A  FILE  OF  EIGENVALUES 
OBTAINED  FROM  THE  SUBROUTINE  LATRV.  THE  USER  IS  REQUIRED  TO  SET 
THE  NUMBER  OF  OBSERVATIONS  AND  THE  NUMBER  OF  INDEPENDENT  VARIABLES. 
ADDITIONALLY,  A  STARTING  FOR  XKRITL  SHOULD  BE  SUPPLIED. 

THE  PROGRAM  WILL  PRINT  VALUES  FOR  THE  INTERGRAL  AND  XKRITL 
AS  THE  PROGRAM  EXECUTES.  IN  THIS  WAY  THE  USER  IS  ABLE  TO 
SEE  THE  HOW  WELL  THE  INTEGRAL  IS  CONVERGING.  IF  THE  INTEGRAL 
IS  CONVERGING  TOO  SLOWLY  THEN  THE  VALUE  OF  XKRITL  SHOULD  BE 
CHANGED 

^AAtt^?0^^a“0Fg%«}0L?®i^100)'L(100)’U(100)'S(100) 

INTEGER  IND.N.iS.K.t.Nft.KP 

open(unit=i£)  , Til£=  mAt 1 1 

OPEN(UNIT=12  FILE='OUT') 

XKRITL=0. 

X K R I T  L=  3 .  40 

ASSIGN  VALUE  TO  ORDER  OF  MATRIX 
N=100 

SET  NUMBER  OF  INDEPENDENT  VARIABLES  (INCLUDE  INTERCEPT  TERM 
K-  6 

IS=N*N 

SET  DESIRED  TRUNCATION  ERROR 
EPS 1=. 0001 

SET  DESIRED  ACCURACY  OF  NUMERICAL  INTEGRATION 
EPS2=. 0001 

SET  DESIRED  ACCURACY  OF  DISTRIBUTION  FUNCTION 
EPS3=. 001 

SET  VALUE  OF  DISTRIBUTION  FUNCTION 

FD=. 05 

SET  INPUT  CODE 

IND=1 

ADJUST  FOR  THE  NUMBER  OF  COLUMNS  OF  THE  MATRICES  'A'  AND 
YX‘  WHICH  ARE  LINEAR  COMBINATIONS  OF  EACH  OTHER 
T=N-K 

GENERATE  THE  'A'  MATRIX 


CALL  SCHMDT 
WRITEflO 
REWIND  lfl 
READ 


,J),  J=1,N),  1=1, N) 


IN 


READ 


’1=1,  IS). 


^OM^U.^E ^TI^E ^IGENVALdES  AND  EIGENVECTORS  OF  MATRIX 


'A' 


'A'  MATRIX  FOR  WHICH  EIGENVALUE  AND  EIGENVECTORS 
.  .WILL  BE  COMPUTED 

.OM0U1 
CALL  LATRV( 

WRITE(12,* 

40  REWIND  12 

READ( 12  *_)(D(  I)  .  1=1, N) 

C  SET  THE  NUMBER  dF  NON-ZERO  EIGENVALUES 

NR=T 

C  SORT  THE  N-K  SMALLEST  EIGENVALUES 

KP=K-1 
DO  10  J=1,T 
L(Jl=D(J+ftP) 

10  CONTINUE 

C  COMPUTE  LOWER  BOUND  FOR  STATISTIC 

CALL  FQUAD(L, NR. FD, EPS1.EPS2 , XKRITL) 

xkritd=abs(xkriTl-^krit6) 

XKRITO=XKRiTL 

IF((ABS(FD-. 05)-EPS3). LE. 0)  GO  TO  30 
IF((FD-.05).  LT.  O)  THEN 
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30 


XKRITL=XKRITL+.  5*XKRITD 
PRINT  *  FD.XKRITL 
GO  TO  46 
END  IF 

IF((FD-.05).GT.  0)  THEN 
XKRITL=XKRITL-.  5*XKRITD 
PRINT  *  FD,XKRITL 
GO  TO  46 
END  IF 

PRINT  *,FD,XKRITL 
STOP 
END 

SUBROUTINE  LATRV 


PURPOSE: 

MATRIX 


FIND  THE  EIGENVALUES  AND  VECTORS  OF  A  REAL  SYMMETRIC 


DESCRIPTION  OF  VARIABLES: 

R=  A  REAL  SYMMETRIC  MATRIX,  DESTROYED  DURING  COMPUTATIONS 
AND  REPLACED  BY  THE  TRANSPOSED  MATRIX  OF  EIGENVALUES 
N=  ORDER  OF  THE  MATRICES 
V=  MATRIX  OF  EIGENVECTORS 
D=  VECTOR  OF  EIGENVALUES  IN  DESCENDING  ORDER 
IND=  INPUT  CODE 

IF  IND=0  EIGENVALUES  AND  EIGENVECTORS  ARE 
COMPUTED 

IF  IND=1  EIGENVALUES  ARE  COMPUTED  ONLY 

METHOD:  DIAGONALIZATION  METHOD  OF  JACOBI 

REFERENCE:  NATIONAL  PHYSICAL  LABORATORIES,  'MODERN  COMPUTING 
METHODS^ 


SUBROUTINE  LATRV  (R,N, V.D, IND) 
DIMENSION  R(  1) ,V(  1 ) ,  D(  1 J 
INDEX(I,J)=I+/6-rrN 

GENERATE  IDENTITY  MATRIX 
IF(IND)  1,1,4 

1  NN=N*N 

IJ=N+1 
DO  2  1=1, NN 

2  V(I)=0 

DO  3  1=1 , NN , I J 

3  V(I)=1 

COMPUTE  FINAL  NORM  EPS 

4  EPS=0 

DO  6  1=1, N 
DO  6  J=1  N 
IJ=INDEX(j  ,1) 

IF(I-J)  5.6.5 

5  EPS=EPS+R( lj)**2 

6  CONTINUE 

I F( EPS-.  0000001)  95,95.7 

7  EPS=SQRT(EPS)*. 0000061/^ 

DREMP=1 

8  L2=0 

DO  80  LR=1,N 
IR=INDEX(0’LR) 

DO  80  LK=Lft,N 
IF(LK-LR)  86,80,10 
10  IK=INDEX(6,LK) 

K1=LR+IK 

K2=LR+IR 

K3=LK+IK 

IF(ABS(ARK)-DREMP)  80,80,20 
20  L2=l 

COMPUTE  SIN  AND  COS 
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26 


C 

C 


;lr,j) 

LK,J) 


40 


45 


50 

80 

90 


95 

100 

110 

115 


120 

130 


135 

140 


ARR= 

AKK 
XU 

XM=ARR-AKK 
SQ=SQRT(XL*XL+XM*XM) 

SIN2T=XL/SQ 
IF(XM)  25.26,26 
SIN2T=-SIN£T  /  N 

SINT2=SINT**2 
C0ST2=C0ST**2 

MODIFY  MATRICES 
DO  40  J=1 ,N 
K4=J+IR 
K5=J+IK 
K6=INDEX 
K7=INDEX 
ARJ=R( K4 
AKJ=R( K5 

R( K4)-AR0*COST+AKJ*SINT 
R(K5)=AKJ*C0ST-ARJ*SINT 
R(K6)=R(K4) 

R(K7)=R(K5) 

CONTINUE 
K7=LK+IR 
R( K1)=0 
K7)=0 

'  K2)=ARR*COST2+AKK*SINT2+ARK*SIN2T 
' K3)=ARR*SINT2+AKK*C0ST2-ARK*SIN2T 
I ND )  45,45,80 
__  50  1=1, N 
K1=I+IR 
K2=I+IK 
ARR=V(K1) 

AKK=V( K2j 

V(K1)=C0$T*ARR+SINT*AKK 
V( K2)-C0ST*AKK-SINT*ARR 
CONTINUE 
IF(L2)  8.90,8 
DREMP=DREfaP*.  1 

COMPARE  THRESHOLD  DREMP  WITH  NORM  EPS 
IF(DREMP-EPSJ  95,8,8 

SORT  EIGENVALUES  AND  EIGENVECTORS 
DO  100  1=1. N 
I J= INDEX ( f , I ) 

D(I )=R( IJ; 

DO  130  1=1, N 
DO  130  J=I , N 
IF(D(J)-D(t)) 

ARR=D(I) 

D( I )=D[J ) 

D(j1=ARR 

IF(IND)  115,115,130 
DO  120  K=1,N 

IJ=INDEXtfc,I) 

IK=INDEX(K  J) 

ARK=V(IJ) 

v(U)=vnK) 

VfIK)=ARK 
CONTINUE 
CONTINUE 

REPLACE  R  BY  TRANSPOSED 
IF(IND)  135,135,150 
DO  140  1=1, N 
DO  140  J=1,N 
IJ=INDEX(J  I) 

IK=INDEX(I  J) 

R(IJ)=V(lk) 


130,130,110 


MATRIX  OF  EIGENVECTORS 


54 


oooooooooooooooooooooo 


150 


RETURN 

END 


SUBROUTINE  FQUAD 


.  LE.  X)  OF 


PURPOSE:  FIND  THE  DISTRIBUTION  FUNCTION  F(X)=P(Q 
QUADRATIC  FORMS  IN  NORMAL  VARIABLES. 

DESCRIPTION  OF  PARAMETERS 

S  =  A  VECTOR  OF  LENGTH  NR  CONTAINING  THE  NON-ZERO 
EIGENVALUES 

NR  =  NUMBER  OF  NON-ZERO  EIGENVALUES 

FD  =  VALUE  OF  THE  DISTRIBUTION  FUNCTION 

EPS1  =  DESIRED  TRUNCATION  ERROR  TU 

EPS2  =  DESIRED  ACCURACY  OF  NUMERICAL  INTEGRATION 

XKRIT  =  VALUE  OF  X 

METHOD:  THE  INTEGRAL  DERIVED  BY  IMHOF  IS  USED  AND  NUMERICALLY 
INTEGRATED  BY  SIMPSONTS  RULE 

REFERENCE:  'COMPUTING  THE  DISTRIBUTION  OF  QUADRATIC  FORMS  IN 
NORMAL  VARIABLES'  BY  J.  P.  IMHOF 
BIOMETRIKA  (1961),  48,  3  AND  4,  P.419 


SUBROUTINE  FQUAD(S, NR, FD, EPS1 , EPS2 , XKRIT) 
DIMENSION  S(I) 

FD=0. 

UB=0. 

C2=0. 

C  1=0. 

FBL=0. 

VINT2=0. 

FD1=0. 

SUMK=0. 

XK=0. 

Si_AM=0. 

DO  6  1=1, NR 
S'" 


XKRIT 


^=fe(5)-> 

c  6  s^M~coi^ufl  u^e^-^iInd^  Vf* Integral  given  the  truncation 

P  FRRHR 

-  •' - -  14472989  +  SLAM)/XK) 


SLAM  +  XK*ALOG(UB) ) 


LT 

0} 

GT. 


0)  GO  TO  20 
GO  TO  9 
0)  GO  TO  9 


10 

11 


14 


15 

16 


UB=EXP(-^AL0G(EPS1*XK)  +  1. 

PRINT  ,UB 

TU=EXP( 1. 14472989  +  ALOG(XK) 

PRINT  *,TU 
VAL=1. ^TU-EPS1 
PRINT  * , VAL 
IF((1.  /'TU-EPS1). 

I Ff (TU-EPS1).  EQ. 

IFCC1-  /TU-EPS1). 

UB=UB  +  5./XK 
GO  TO  7 

COMPUTATION  OF  THE  INTEGRAND  GIVEN  THE  VALUE  OF  U 
I F(U)  15,15,11 
TETA=Q. 

RH0=1. 

DO  14  1=1, NR 
C1=S(I)*U 

C2=l.  +  Cl**2  /  N 

TETA=TETA  +  ATAN(Cl) 

RH0=RH0*(C2**.  25) 

TETA=. 5^(TETA) 

FBL=SIN(TETA)/(RHO*U) 

GO  TO  18 
FBL=0. 

DO  16  1=1. NR 
FBL=FBL+Su) 

FBL=.  5*( FBL-XKRIT) 
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18  GO  TO  (21.22,23.24,25).  KSKIP 

EVALUAYkJn  OF  fHE  TNTEGRAL  BY  SIMPSON'S  RULE 

20  RANGE=UB 

MH=1 

U=RANGE*.  5 
KSKIP=1 
GO  TO  10 

21  SUMK=FBL*RANGE*2./3. 

U=0. 

KSKIP=2 
GO  TO  10 

22  V I NT2=SUMK+FBL* RANGE/6. 

U=UB 
KSKIP=3 
GO  TO  10 

23  VINT2=VINT2  +  FBL*RANGE/6. 

FD=. 5-. 318309886*VINT2 
DO  28  NIT=1 , 14 
FD1=FD 

VINT2=(VINT2-SUMK*.  5)*.  5 

MH=2*MH 

STEP=RANGE/MH 

U=STEP*.  5 

KSKIP=4 

GO  TO  10 

24  SUMK=FBL 

IF(MH.  EQ.  2)  THEN 

U=U+STEP 

KSKIP=5 

GO  TO  10 

END  IF 

DO  25  K=2 ,MH 
U=U+STEP 
KSKIP=5 
GO  TO  10 
SUMK=SUMK+FBL 
SUMK=SUMK*STEP*2./3. 

V I NT  2=V I NT 2  +  SUMK 
FD=.  5  -  ,318309886*VINT2 
IF(NIT-3)  28,28.27 
IF(ABS(FDl-FD)-Ef>S2)  29,28,28 
CONTINUE 
PAUSE  1234 
RETURN 
END 

SUBROUTINE  SCHMDT 


25 


27 

28 

29 


PURPOSE:  GENERATE  THE 
COMPUTED 


'A'  MATRIX  FOR  WHICH  EIGENVALUES  ARE 


DESCRIPTION  OF  PARAMETERS 

A  =  MATRIX  FOR  WHICH  EIGENVALUES  ARE  COMPUTED 
N  =  ORDER  OF  MATRIX 


SUBROUTINE  SCHMDT(A,N) 
REAL  A( 100 , 100) 

DO  10  1=1, N 
DO  20  J=1  N 
A(I.J)=0. 

20  CONTINUE 
10  CONTINUE 
A  1,1)=2. 

Af  1  2)=-l. 

A(  1  5  =-l. 

DO  3o  1=2 , N 
J=I 

IFfI.LE.4)  THEN 
A(I,J)=3. 
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;j+n)=-i. 


IF((  I.  GE.  (N-3) ).  AND.  (I.  LT.  N))  THEN 
A(  I  ,J)=3. 
a  i;ro+i)  =-i. 


IF(I.EQ.N)  THEN 
A(I,J}=2. 

Afl  (J-l))=-l. 

A(I  (J-4))=-l. 

END  IF 

IF((  I.  GT.  4).  AND.  (I.  LT.  (N-3)))  THEN 

J-JH 


CONTINUE 

RETURN 

END 
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THIS  PROGRAM  COMPUTES  THE  UPPER  BOUND  FOR  THE  SCHMIDT  STATISTIC 


LE  MAT  IS  A  FILE  OF  EIGENVALUES 
RV.  THE  USER  IS  REQUIRED  TO  SET 
HE  NUMBER  OF  INDEPENDENT  VARIABLES. 


C 

C 


C 

C 


OPEN( _  _ 
OPENfUNIl 
XKRIT  ‘ 


USING  IMHOFF'S  INTEGRAL  THE  F 
OBTAINED  FROM  THE  SUBROUTINE  LA 

THE  NUMBER  OF  OBSERVATIONS  AND  ....  _  _ 

ADDITIONALLY,  A  STARTING  FOR  XKRITU  SHOULD  BE  SUPPLIED. 

THE  PROGRAM  WILL  PRINT  VALUES  FOR  THE  INTERGRAL  AND  XKRITU 
AS  THE  PROGRAM  EXECUTES.  IN  THIS  WAY  THE  USER  IS  ABLE  TO 
SEE  THE  HOW  WELL  THE  INTEGRAL  IS  CONVERGING.  IF  THE  INTEGRAL 
IS  CONVERGING  TOO  SLOWLY  THEN  THE  VALUE  OF  XKRITU  SHOULD  BE 
CHANGED 

BSMts:aiiisiot)'“(ioou(ito)'"(M>'s(ioi) 

INTEGER  IND.N.lS.K.t.Nft.KP 
.unit=iMil£='mAt') 

UNIT=12  FILE='OUT') 

ru=o. 

XKRITU=3  80 

ASSiGN  VALUE  TO  ORDER  OF  MATRIX 

N=20 

SET  NUMBER  OF  INDEPENDENT  VARIABLES  (INCLUDE  INTERCEPT  TERM 
K=5 

j  ^ ^  N 

SET  DESIRED  TRUNCATION  ERROR 
EPS1=. 0001 

SET  DESIRED  ACCURACY  OF  NUMERICAL  INTEGRATION 
EPS2=  0001 

SET  DESIRED  ACCURACY  OF  DISTRIBUTION  FUNCTION 
EPS3=. 001 

SET  VALUE  OF  DISTRIBUTION  FUNCTION 

FD=. 05 

SET  INPUT  CODE 

IND=1 

ADJUST  FOR  THE  NUMBER  OF  COLUMNS  OF  THE  MATRICES  'A'  AND 
YX'  WHICH  ARE  LINEAR  COMBINATIONS  OF  EACH  OTHER 
T=N-K 

GENERATE  THE  'A'  MATRIX 


CALL  SCHMDT 
WRITE! 10  *' 
REWIND  1(5 
READ  IN 


,J),  J=1,N),  1=1, N) 


■A1  MATRIX  FOR  WHICH  EIGENVALUE  AND  EIGENVECTORS 
WILL  BE  COMPUTED 
READ  (10  *5  (R(I).  1=1, IS) 

OMI^UTE  THE  EIGENVALUES  AND  EIGENVECTORS  OF  MATRIX  'A' 
CALL  LATRVfR.N.V.D.IND) 


wrTt  e^12  ,*)(i5(  i )  1^1=1^) 

40  REWIND  12 

READ( 12  *_)(D( I) .  1=1. N) 

C  SET  THE  NUMBER  dF  NON-ZERO  EIGENVALUES 

NR=T 

C  SORT  THE  N-K  LARGEST  EIGENVALUES 

DO  20  1=1, T 

um=Dm 

20  CONTINUE 

C  COMPUTE  LOWER  BOUND  FOR  STATISTIC 

CALL  FQUADTu. NR, FD.EPS1.EPS2. XKRITU) 

C  SUCCESSIVELY  HAEf  XKRITU  fO  FIND  THE  5  PERCENT 

C  ,  SIGNIFICANCE  POINT 

XKRITD=ABS(XKRITU-XKRITO) 

XKRITO=XKRITU 

IF(!ABS(FD-.  05)-EPS3).  LE. 0)  GO  TO  30 
IF(fFD-.  05).  LT.  O)  THEN 
XKRITU=XKRITU+. 5*XKRITD 
PRINT  *.FD, XKRITU 
GO  TO  46 
END  IF 

IF((FD-.05).GT.  0)  THEN 
XKRITU=XKRITU-. 5*XKRITD 
PRINT  *AFD, XKRITU 
GO  TO  46 
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END  IF 

30  PRINT  *,FD,XKRITU 
STOP 
END 

SUBROUTINE  LATRV 


PURPOSE: 

MATRIX 


FIND  THE  EIGENVALUES  AND  VECTORS  OF  A  REAL  SYMMETRIC 


DESCRIPTION  OF  VARIABLES: 

R=  A  REAL  SYMMETRIC  MATRIX,  DESTROYED  DURING  COMPUTATIONS 
AND  REPLACED  BY  THE  TRANSPOSED  MATRIX  OF  EIGENVALUES 
N=  ORDER  OF  THE  MATRICES 
V=  MATRIX  OF  EIGENVECTORS 
D=  VECTOR  OF  EIGENVALUES  IN  DESCENDING  ORDER 
IND=  INPUT  CODE 

IF  IND=0  EIGENVALUES  AND  EIGENVECTORS  ARE 
COMPUTED 

IF  IND=1  EIGENVALUES  ARE  COMPUTED  ONLY 

METHOD:  DIAGONALIZATION  METHOD  OF  JACOBI 

REFERENCE:  NATIONAL  PHYSICAL  LABORATORIES,  'MODERN  COMPUTING 
METHODS' 


SUBROUTINE  LATRV  (R,N,V.D,IND) 
DIMENSION  R(1 ) . V( 1)  ,D( l) 
INDEXa.J)=I+fJ-lJ*N 

GENERATE  IDENTITY  MATRIX 
IF(IND)  1,1,4 

1  NN=N*N 

IJ=N+1 
DO  2  1=1, NN 

2  V(I)=0 

DO  3  1=1 ,NN , IJ 

3  V(I)=1 

COMPUTE  FINAL  NORM  EPS 

4  EPS=0 

DO  6  1=1, N 
DO  6  J=1,N 
IJ=INDEX(J,I) 

IF(I-J)  5,6,5 

5  EPS=EPS+R(lj)**2 

6  CONTINUE 

IF( EPS-.  0000001)  95,95.7 

7  EPS-SQRT(EPS)*.  00000t)l/N 

DREMP=1 

8  L2=0 

DO  80  LR=1,N 

ir=index(o!lr) 

DO  80  LK=LAaN 
IF(LK-LR)  86,80,10 
10  IK=INDEX(6,LK) 

K1=LR+IK 

K2=LR+IR 

K3=LK+IK 

IF(ABS(ARK)-DREMP)  80,80,20 
20  L2=l 

COMPUTE  SIN  AND  COS 
ARR=R(K2) 

AKK=R(K3) 

XL=2*ARK 
XM=ARR-AKK 
SQ=SQRT(XL*XL+XM*XM) 
SIN2T=xL/SQ 
IFCXM)  25,26,26 
25  SIN2T=-SIN^T 
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;lr,j) 

LK,J) 


«RNT^^sxW/SQ)/2) 

SINT2=SINT**2 
COST2=COST**2 

C  MODIFY  MATRICES 

DO  40  J=1,N 
K4=J+IR 
K5=J+IK 
K6=INDEX( 

K7= I NDEX ( 

ARJ=R(K4 
AKJ=R(K5l 

R( K4)-ARJ*C0ST+AKJ*SINT 
R(K5)=AKJ*C0ST-ARJ*SINT 
R( K6)=R( K4) 

R(K7)=R(K5) 

40  CONTINUE 
K7=LK+IR 
R(K1)=0 
R(K7)=0 

R( K2)=ARR*C0ST2+AKK*SINT2+ARK*SIN2T 
RLK3)=ARR*SINT2+AKK*COST2-ARK*SIN2T 
IF(IND)  45,45,80 
45  DO  50  1=1, N 
K 1= I + I R 
K2=I+IK 
ARR=V(K1) 

AKK=V(K2) 

V(K1)=C0ST*ARR+SINT*AKK 
50  V(K2l=C0ST*AKK-SINT*ARR 
80  CONTINUE 

IF( L21  8,90,8 
90  DREMP=DREfoP*!  1 

C  COMPARE  THRESHOLD  DREMP  WITH  NORM  EPS 

IF(DREMP-EPS)  95,8,8 

C  SORT  EIGENVALUES  AND  EIGENVECTORS 

95  DO  100  1=1. N 
IJ=INDEXn,I) 

100  D(I)=R(IJ) 

DO  130  1=1, N 
DO  130  J=I ,N 

IF(D(J)-D(i))  130,130,110 
110  ARR=D( I) 

D(I)=DfJ) 

D(J)=ARR 

IF(IND)  115,115,130 
115  DO  120  K=1 ,N 
IJ=INDEX(K, I) 

IK=INDEX(K  J) 

ARK=V( IJ) 

V(IJ)=V[IK) 

v(ik1=ark 

120  CONTINUE 
130  CONTINUE 

C  REPLACE  R  BY  TRANSPOSED 

IF(IND)  135,135,150 
135  DO  140  1=1. N 
DO  140  J=1,N 
IJ=INDEX(J  I) 

IK=INDEX(I,J) 

140  R(IJ)=V(IK) 

150  RETURN 
END 

SUBROUTINE  FQUAD 

PURPOSE:  FIND  THE  DISTRIBUTION  FUNCTION  F(X)=P(Q  . LE.  X)  OF 
QUADRATIC  FORMS  IN  NORMAL  VARIABLES. 

DESCRIPTION  OF  PARAMETERS 


MATRIX  OF  EIGENVECTORS 
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S  =  A  VECTOR  OF  LENGTH  NR  CONTAINING  THE  NON-ZERO 
EIGENVALUES 

NR  =  NUMBER  OF  NON-ZERO  EIGENVALUES 

FD  =  VALUE  OF  THE  DISTRIBUTION  FUNCTION 

EPS1  =  DESIRED  TRUNCATION  ERROR  TU 

EPS2  =  DESIRED  ACCURACY  OF  NUMERICAL  INTEGRATION 

XKRIT  =  VALUE  OF  X 

METHOD:  THE  INTEGRAL  DERIVED  BY  IMHOF  IS  USED  AND  NUMERICALLY 
INTEGRATED  BY  SIMPSON'S  RULE 

REFERENCE:  'COMPUTING  THE  DISTRIBUTION  OF  QUADRATIC  FORMS  IN 
NORMAL  VARIABLES'  BY  J. P.  IMHOF 
BIOMETRIKA  (1961),  48,  3  AND  4,  P.419 


SUBROUTINE  FQUAD(S , NR , FD, EPS1 , EPS2 , XKRIT) 

DIMENSION  S(l) 

FD=0. 

UB=0. 

C2=0. 

C1=0. 

FBL=0. 

VINT  2=0 . 

FD1=0. 

SUMK=0. 

XK=0. 

SLAM=0. 

DO  6  1=1, NR 

XK=Ik+( 5  _XKRIT 

6  SLAM=SLAM+.  5*ALOG(ABS(S( I))) 

C  COMPUTE  UPPER-BOUND  OF  INTEGRAL  GIVEN  THE  TRUNCATION 

C  ERROR 

UB=EXP(-(AL0G(EPS1*XK)  +  1.14472989  +  SLAM)/XK) 

C  PRINT  *  UB 

7  TU=EXP(’l.  14472989  +  ALOG(XK)  +  SLAM  +  XK*ALOG(UB)) 

C  PRINT  , TU 

VAL=1.>TU-EPS1 
C  PRINT  *,VAL 

IF((1.  ^TU-EPS1). LT.  0)  GO  TO  20 
I Fu  TU-EPS1).  EQ.  0}  GO  TO  9 
IFCCl- /TU-EPSl).  GT.  0)  GO  TO  9 
9  UB=UB  +  5./XK 

GO  TO  7 

COMPUTATION  OF  THE 
IF(U)  15,15,11 
TETA=0. 

RHO=l. 

DO  14  1=1, NR 
C1=S(I)*U 

C2=l.  +  C 1**2  ,  v 
TETA=TETA  +  ATAN(Cl) 

RHO=RHO*(C2**.  25) 

TETA=.  5*fTETA) 

FBL=SIN(TETA)/(RHO*U) 

GO  TO  18 
FBL=0. 

DO  16  1=1. NR 
FBL=FBL+S(1) 

FBL=. 5*(FBL-XKRIT)  % 

G°  TE°V AL2U A’T I CfN  3Q F  V H E5  f N T ESG R'AL  BY  SIMPSON'S  RULE 
RANGE=UB 
MH=1 

U=RANGE*. 5 
KSKI P=1 
GO  TO  10 

21  SUMK=FBL*RANGE*2.  /3. 


10 

11 


14 


15 

16 
18 
20 


INTEGRAND  GIVEN  THE  VALUE  OF  U 
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22 


23 


U=0. 

KSKIP=2 
GO  TO  10 

VINT2=SUMK+FBL*RANGE/6. 
U=UB 
KSKIP=3 
GO  TO  10 

VINT2=VINT2  +  FBL*RANGE/6. 
FD=.  5-.  318309886*VINT2 
DO  28  N I T= 1 , 14 
FD1=FD 

VINT2=(VINT2-SUMK*.  5)*.  5 

MH=2*MH 

STEP=RANGE/MH 


24 


U=S 

KSK 

GO 


THEN 


25 


27 

28 

29 


TEP*. 
tP=4 

.  ro  10 

SUMK=FBL 
IF(MH.  EQ.  2) 

U=U+STEP 
KSKIP=5 
GO  TO  10 
END  IF 

DO  25  K=2 , MH 
U=U+STEP 
KSKI P=5 
GO  TO  10 
SUMK=SUMK+FBL 
SUMK=SUMK*STEP*2.  /3. 
VINT2=VINT2  +  SUMK 
FD=.  5  -  .  318309886*VINT2 
I FC NIT— 3^  28,28,27 
IF(ABS(FDl-r 
CONTINUE 
PAUSE  1234 
RETURN 
END 

SUBROUTINE  SCHMDT 


CO .CO .c/  v 

-FD)-Ef>S2)  29,28,28 


PURPOSE:  GENERATE  THE 
COMPUTED 


•A’  MATRIX  FOR  WHICH  EIGENVALUES  ARE 


DESCRIPTION  OF  PARAMETERS 

A  =  MATRIX  FOR  WHICH  EIGENVALUES  ARE  COMPUTED 
N  =  ORDER  OF  MATRIX 


SUBROUTINE  SCHMDT(A,N) 
REAL  A(IOO.IOO) 

DO  10  1=1, N 
DO  20  J=1 ,N 
A(I.J)=0. 

20  CONTINUE 
10  CONTINUE 
A( 1,1 )=2. 

Afl  2)=-l. 

Afl’5)=-1. 

DO  50  I=2,N 

IFfl.LE.  4)  THEN 

=-ll 
=-l. 


AND.  (I.  LT.N))  THEN 
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30 


Stf-K-4”-1- 

IF(I.EQ.N)  THEN 


IF(( I.  GT.  4).  AND.  (I.  LT.  (N-3)))  THEN 


RETURN 
END 
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Program  Stats  is  an  APL  simulation  which  was  used  to  generate  data  for  the 
comparison  of  the  estimates  of  and  the  ability  of  the  Schmidt  statistic  as  part  of  the 
AR(1,4)  procedure  and  the  Durbin- Watson  statistic  as  part  of  the  two-step  procedure 
to  detect  the  AR(1,4)  disturbance. 

V  STATS  ;ADliAD;AWlBi  AW1 1 AW2  ;  AW  3  ;  AWnB ;  AW*+  ;  AW5  ;  AW ;  AS ;  PW 
[  1  ]  n  ENTER  VALUES  FOR  RHO 1  AND  RHOn 
[  2  ]  '  ENTER  A  VALUE  FOR  RHO  1 ' 

[3]  RHOl+O 

[  4  ]  '  ENTER  A  VALUE  FOR  RHO 4  ' 

[5]  RHOV+Q 

[  6  ]  '  ENTER  NUMBER  OF  OBSERVATIONS  • 

[7]  NOBS* □ 

[  8  ]  '  ENTER  VARIANCE  OF  NORMAL  DISTRIBUTION ' 

[9]  SIG2+U 

[10]  A  INITIALIZE  BOUND  FOR  STATISTICAL  TESTS 

[11]  DWL+l. 65 

[12]  DWU+l . 69 

[13]  DWAL*-1 .  59 

[14]  DWAU+1 . 63 

[15]  DSMTL+3 .45 

[16]  DSMTU+3. 53 

[  1 7]  fi  INITIALIZE  OUTER  LOOP  VECTORS  AND  COUNTERS 

[18]  RHOAV l-e-40p0 

[19]  i?ff0AP4S«-4OpO 

[20]  RHOAV^+HOpO 

[21]  Z?l«-40p 0 

[22]  i?4S<-40 p 0 

[23]  £4-e-40p0 

[24]  DWCNT<rO 

[25]  DSTCNT<rO 

[26]  DWACNT+ 0 

[27]  DWINCT^O 
[2  8]  STINCT<rO 

[29]  DWAICT+Q 

[30]  n  GENERATE  DURBIN -WATSON  'A'  MATRIX 
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[31]  NO+NOBS 

[32]  ADlB*-(  (  (NO -2  )p0),”l,2,~l)),(  {NO -2  )pO  )  ,  ”l ,  1 

[33]  AD1<-1,~1,  ((  (NO*  (.NO-2)  )+  (NO-2)  )pADlB 

[34]  AD+(NO ,NO)pADl 

[3  5]  o  GENERATE  THE  WALLIS  'A'  MATRIX 
[36]  AWlB*-(  (  (NO -5  )-3  )pO  ) 

[3  7]  AW  1C*-  (  (  (  (NO -5  )  +  l)p0),l,0,0,0,~l))  ,AW1B  ) 

[3  8]  AJ71-«-l  ,0,0,0,-l,((((3  *N0  )+3  )  pAWlC 

[3  9]  AW2*-(  (  (NO -8  )*N0  )  —  9  )  p  (  (  (NO -8  )p0),_l,0,0,0,2,0,0,0,~l) 

[40]  AW3*-~  1 , 0,0, 0,2,0, 0,0, ~1  ,AW2 

[41]  AWUB+  (((3xNO  )  +  3  )  )p  ((.((NO- 5)+l)p0),  "1,0, 0,0,1) 

[42]  AWU*-(  (  (NO- 5  )-3  )p0  )  ,  ~1 ,  0 , 0 , 0 , 1  ,AWVB 

[43]  AW5*-AW1 ,AW3 ,  AJ74 

[44]  AW*-(NO  ,NO)pAW5 

[45]  fl  GENERATE  THE  SCHMIDT  'A'  MATRIX 

[46]  AS+AD+AW 

[47]  A  OUTE^LOOP  FOR  REPLICATION 

[48]  COUNT 0+0 

[49]  OUTER :COUNTO+COUNTO+ 1 

[5  0]  fl  INITIALIZE  INNER  LOOP  VECTORS 

[51]  RHOHAT1*-10 p  0 

[52]  RHOHATnS*- lOpO 

[53]  RHOHATn+lOpO 

[  54]  fl  INNER  LOOP  TO  CALCULATE  RHO 1 ,  RHOn  AND  DETERMINE  DETECTION 
[5  5]  COUNT  I <-0 

[56]  INNER: COUNTI+COUNTI+l 

[57]  DETDW+0 

[58]  DETSMT+o 

[59]  DETWA+0 

[60]  INSTCT*-0 

[61]  INDUCT*- 0 

[62]  INDWAT+0 

[63]  fl  GENERATE  THE  RANDOM  ERROR 

[64]  V*-NOBS  NORRAND (0 , SI G2) 

[65]  E+NOBSpO 
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[66]  JV+4 

[67]  A  CALCULATE  CORRECTION  TERMS 

[68]  c7>(l-(fiffO4*2))*0.5 

[69]  I+-(RH01*J) 

[70]  G+RHOUxI 

[71]  F+(1+(RH01*2  )-(RROn*2  )-  (J* 2  )  )*0 . 5 

[72]  El<r-(RH01*F) 

[73]  D<—  (I*G)*F 

[74]  C+(1  +  (.RH01*2  )-  (RHOU*2  )-  (El*2  )  )*0. 5 

[75]  Bl<-(i?ff01-(0xEl))  +  C 

[76]  A<r(l-(RHOU*2  )-  (C*2  )-(£>*2  )-(Bl*2  )  )*0 . 5 
[7  7]  A  CALCULATE  THE  FIRST  FOUR  ERROR  TERMS 

[78]  £[1]«-7[1]*4 

[79]  £[2]<-(7[2]-(Blx£[l]  ))*C 

[80]  £[3]«-(7[3]-(Bx£[l]  )-(£lx£[2]  ))*F 

[81]  E[4]«-(y[4]-(GxE[l]  )-(Jxff[3]  ))*J 

[82]  MM:N+N+ 1 

[83]  E[JV]«-(flff01x£[iV-l]  )+(B«04xF[fl/-4]  )+7[A/] 

[84]  +(N<NOBS)/MM 

[85]  £<-(W0BS,l)p£ 

[86]  r  CBJVBBArB  THE  INDEPENDENT  VARIABLES 

[87]  Xl«-OV0BS,l)pOV0BSpl) 

[88]  X2+(NOBS ,  1  )p  (.NOBS  NORRAND  0  0.0625  ) 

[89]  X«-X1,X2 

[  9  0  ]  R  GENERATE  THE  TRUE  BETA 
[91]  BETAT+  2  1  p  1  1 

[9  2]  fl  GENERATE  THE  DEPENDENT  VARIABLES 

[93]  Y+(X+.xBETAT)+E 

[94]  fl  LEAST  SQUARES  ESTIMATE  OF  BETA 

[95]  B+Y1X 

[96]  fl  GENERATE  THE  RESIDUALS 

[97]  EMW-(X+.xB) 

[98]  fl  CALCULATE  THE  DURBIN -WATSON  STATISTIC 

[99]  DW<- ,  (  (§EHAT)  + .  x  (AD+.  *EHAT)  )*  (  (§EHAT)  + .  *EHAT) 
[10  0]  fl  CALCULATE  THE  SCHMIDT  STATISTIC 
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[101]  DS+ ,  ( ( §EHAT  )  +  .  X  QS+  .  *EHAT )  )  *  ( ( §EHAT )  +  .  x£ff  Ar ) 

[102]  n  COMPARE  TABLED  VALUE  FOR  DURBIN -WATSON  STATISTIC  TO 
[10  3]  fi  THE  ABOVE  CALCULATED  VALUE  FOR  THE  DURBIN -WAT SON 

[104]  DETDW+iDWZDWL) 

[105]  DWCNT+DWCNT+DETDW 

[106]  INDWCT+ ( (DWL<DW )a (DWU>DW ) ) 

[107]  DWINCT+DWINCT+INDWCT 

[10  8]  fi  COMPARE  TABLED  VALUE  FOR  SCHMIDT  STATISTIC 

[109]  A  TO  ABOVE  CALCULATED  VALUE  FOR  SCHMIDT  STATISTIC 

[110]  DETSMT+ ( DSZDSMTL ) 

[111]  DSTCNT+DSTCNT+DETSMT 

[112]  INSTCT+((DSMTL<DS')k(.DSMTU>DS)) 

[113]  STINCT+STINCT+INSTCT 

[114]  R  CALCULATE  THE  DURBIN -WATSON  ESTIMATE  FOR  RHOl 

[115]  RHOHAT11COUNTI1+1- (DW*2) 

[116]  fi  CALCULATE  THE  SCHMIDT  ESTIMATE  FOR  RHO^ 

[117]  RHVS1+2-RHOHAT11COUNTI1 

[118]  RHnS2+DS*2 

[119]  RHOHAT^S [COUNT I] +RHHS1 -RHnS2 

[120]  fi  GENERATE  TRANSFORM  MATRIX  FOR  FIRST  ORDER  PROCESS 

[121]  PD+,  (l-(RHOHATlLCOUNTn*2)')*0.5 

[122]  PW+( (NO*NO)-l )p ( ( (NO- 1 )p0 ) , (-RHOHAT1 [ COUNTI1 ) ,1 ) 

[123]  PDW+(NO ,NO)pPD ,PW 

[124]  fi  TRANSFORM  RESIDUALS 

[125]  ESTAR+PDW+.xEHAT 

[126]  fi  CALCULATE  THE  WALLIS  STATISTIC 

[127]  DWA+ , ( ( §ESTAR ) + . x ( AW+ . *ESTAR ) ) * ( ( SESTAR ) + . xESTAR ) 
[12  8]  fl  COMPARE  TABLED  VALUE  OF  WALLIS  STATISTIC 

[129]  fl  TO  CALCULATED  VALUE  OF  WALLIS  STATISTIC 

[130]  DETWA+(DWAZDWAL) 

[131]  DWACNT+DWACNT+DETWA 

[13  2]  INDWAT<r (  (DWAL<DWA)a(DWAU>DWA)) 

[133]  DWAICT+DWAICT+INDWAT 

[134]  a  CALCULATE  THE  WALLIS  ESTIMATE  FOR  RHOU 

[135]  RHOHATmcOUNTIl+1- (DWA*2) 
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[136]  +(C0UNTI<5 )  /INNER 

[137]  REOAVllCOUNTOl+(+/RROEATl )*5 
[13  8]  RHOA  V  45  [  CO  UN  TO  ]  <-  (  +  /RHOHA  T^S)*5 
[13  9]  RE  OA  74  ICO  UNTO  ]-«-(  +  /RE  OR  A  T4  )  t  5 

[140]  +(COUNTO<nO ) /OUTER 

[141]  CNTDW+DWCNT 

[142]  CNTDST+DSTCNT 

[143]  CNTDWA +DWA CNT 

[144]  R1+RHOAV1 
[14  5]  RHS+REOAV^S 
[146]  RU+REOAVn 

V 
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Program  Regress  is  an  APL  simulation  which  generates  data  for  the  comparison 
of  the  regression  parameters  from  both  the  two-step  and  AR(1,4)  procedures. 

V  REGRESS  -,AD  iADUAWIB  lAWl  iAW2  ;AW3  iAWUB  lAWn  ;AW5  iAW iAS 
[  1  ]  A  ENTER  VALUES  FOR  RHO 1  AND  RHOn 
[2  ]  '  ENTER  A  VALUE  FOR  RHOl ' 

C  3 ]  RHOl+U 

[4]  '  ENTER  A  VALUE  FOR  RHOU  ' 

[5]  RHOU+U 

[  6  ]  '  ENTER  NUMBER  OF  OBSERVATIONS ' 

[7]  NOBS+U 

[  8  ]  1  ENTER  VARIANCE  OF  NORMAL  DISTRIBUTION ' 

[9]  SIG2+U 

[10]  DSMTL+3.H 5 

[11]  DSMTU+ 3.53 

[12]  A  GENERATE  DURBIN -WATSON  '  A  '  MATRIX 

[13]  NO+NOBS 

[14]  ADlB+( ( (NO -2 )p0),~l,2,~l)), ( (NO-2 )p 0 ) , "l , 1 

[15]  AD1+1,~1,  (((NO*(NO-2))+(NO-2))pADlB 

[16]  AD+(NO ,NO)pADl 

[17]  A  GENERATE  THE  WALLIS  '  A  '  MATRIX 

[18]  AWlB+(((NO-5)-3)pO) 

[19]  AWlC+(  (  (  (.NO -5  )+l)p0),l,0,0,0,”l))  ,AW1B  ) 

[20]  ^J71-«-l  ,0,0,0,”l,(((  (3*N0  )+3  )pAWlC 

[21]  AW  2+  (  (  (NO -3  )*NO )  —  9  )  p  (  (  (NO -8  )p0),~l,0,0,0,2,0,0,0,~l) 

[22]  AW3+~1, 0,0, 0,2, 0,0,0,  ~1,AW2 

[23]  AWUB+( ((3*N0 )+3 ) )p ( ( ( (W0-5 )+l )p0 ) , "l , 0,0, 0,1) 

[24]  AWn<r(((NO-5)-3  )p0  >,"1,0,0,0,1,711/45 

[25]  AW5+AW1 ,AW3 ,AWU 

[26]  AW+(NO ,NO)pAW5 

[27]  A  GENERATE  THE  SCHMIDT  'A'  MATRIX 

[28]  AS+AD+AW 

[  2  9  ]  A  GENERATE  TRANSFORM  MATRIX  FOR  FIRST  ORDER  PROCESS 

[30]  NO+NOBS 

[31]  PD+,(1~(RH01*2))*0 .5 

[32]  PW+((NO*NO)-l)p(((NO-l)pO) ,(~RHOl) ,1) 
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[33]  PDW+(NO ,NO)pPD ,PW 

[34]  n  GENERATE  TRANSFORM  FOR  FOURTH  ORDER  PROCESS 

[35]  Pl«-(l-(tf//04*2))*0.5 

[36]  P41«-P1,  ((/VO-l)pO) 

[37]  P42*-0,P1  ,  (  (NO-2  )pO  ) 

[38]  P43«-0 , 0 ,P1 , (  (NO -3  )p0  ) 

[39]  P44«-0,0,0,P1,  ((/V0"4)p0) 

[40]  P4  5«-  (  (/VO-4  )x/V0)p((-BB04),  0,0, 0,1,  ((/VO-4  )p0)) 

[41]  P4«-(JV0  ,/VO  ) p  (P41  ,P42  ,P4  3  ,P44  ,P4  5  ) 

[42]  n  INITIALIZE  MATRICES  AND  VECTORS  FOR  OUTER  LOOP 

[43]  BLSIAV+VOpO 

[44]  BLSSAV+^OpO 

[45]  B2SIAV+U0p0 

[46]  B2SS71V«-40p0 

[47]  B14J717«-40p0 

[48]  B14S717«-40p0 

[49]  DSBCNT+0 

[50]  DSACNT+ 0 

[51]  STIBCT+ 0 

[52]  STANCT+ 0 

[  53  ]  A  OUTER  LOOP  FOR  REPLICATION 

[54]  COUNTO<- 0 

[55]  OUTER : COUNTO+COUNTO + 1 

[  56]  p  INITIALIZE  MATRICES  AND  VECTORS  FOR  INNER  LOOP 

[57]  BLSI+ 5p0 

[58]  BLSS+ 5p0 

[59]  B2Si>5pO 

[60]  B2SS+5pO 

[61]  B14J«-5p0 

[62]  P140^-5p 0 

[63]  P  INNER  LOOP  TO  CALCULATE  RHO 1  ,  RHOn  AND  DETERMINE  DETECTION 

[64]  COUNT 1+ 0 

[65]  INNER : COUNTI+COUNTI+1 

[66]  P  GENERATE  THE  RANDOM  ERROR 

[67]  V+NOBS  NORRANDiO  ,SIG2) 
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[68]  E+NOBS pO 

[69]  N+ 4 

[70]  o  CALCULATE  CORRECTION  TERMS 

[71]  -  (RHOU*2  )  )*0 . 5 

[72]  I-e--(Ptf01*c7) 

[73]  C+RHOUxI 

[74]  F«-(l  +  (m?l*2  )-  (FF04*2  )-  (1*2  ) )*0 . 5 

[75]  E1+-(RH01*F') 

[76]  D+-U*G)*F 

[7  7]  C-e-(l  +  (FF01*2)-(FP04*2)-(Fl*2))*O.5 

[78]  Bl*-(RH01-(DxEl))±C 

[79]  71-Kl  -  (FF04*2  )-(G*2 )-  (D* 2  )-(B  1*2  )  )*0 . 5 

[  8  0  ]  R  DETERMINE  TRANSFORM  MATRIX  FOR  1 , 4  PROCESS 

[81]  P141-e-d,  ((iVO-l)pO) 

[82]  P14  2<-B1  ,C,(  (NO -2  )p0  ) 

[83]  P143*P,F1,F,  (  (iVO-3  )p0  ) 

[84]  P144-«-(7, 0 ,1 ,J ,  (  (iVO-4  )p0  ) 

[8  5]  P14  5<-((W-4)x^O)p((-PP04),0,0,  ( -FP01  )  ,  1 ,  (  (iV0-4  )p  0  )  ) 

[86]  PlU+iNO ,NO)p (PI 41 ,P142,P143 ,P144,P145 ) 

[87]  R  CALCULATE  THE  FIRST  FOUR  ERROR  TERMS 

[88]  Elll<-Vlll*A 

[89]  F[2]«-(7[2]-(BlxF[l]  ))*C 

[90]  F[3]-e-(7[3]  ~(Pxp[i]  )-(Flx£[2]  ))*F 

[91]  F[4]-e-(7[4]  ~(CxF[l]  )-(JxF[3]  ))*<7 

[92]  MM:N+N+1 

[93]  EIN1+(RH01*ELN-1]  )  +  (RHOU*ElN-Ul  )+7[iV] 

[94]  ->(N<NOBS)/MM 

[95]  E+(NOBS,l )pF 

[96]  p  GENERATE  THE  INDEPENDENT  VARIABLES 

[97]  XI-*- (.NOBS ,  1  )p  (NOBSpl  ) 

[98]  X2-diV0FS,l)p(W0flSWFFAWP  0  0.0625) 

[99]  X+X1,X2 

[100]  R  GENERATE  THE  TRUE  BETA 

[101]  BET  AT-*-  2  1  p  1  1 

[10  2]  R  GENERATE  THE  DEPENDENT  VARIABLES 
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[103]  Y*(.X+.xBETAT)+E 

[104]  q  LEAST  SQUARES  ESTIMATE  OF  BETA 

[105]  BLS+YEX 

[106]  BLSIZCOUNTI1+BLSZ 1;] 

[107]  BLSSZCOUNTIl+BLSl 2;] 

[10  8]  n  rffE  i?ES  JZW71BS 

[109]  £Mr-e-J-(X+.xSLS) 

[110]  r  USE  FIRST  ORDER  P  MATRIX  TO  TRANSFORM  Y  ,X ,  AND  EH  AT 

[111]  r TRANSFORM Y 

[112]  YSTAR1+PDW+ . xj 

[113]  o  TRANSFORM  X 

[114]  XSMfllfPZW+.xX 

[115]  r  TRANSFORM  RESIDUALS 

[116]  ESTAR1+PDW+. xEHAT 

[117]  o  USE  FOURTH  ORDER  P  MATRIX  TO  TRANSFORM  YSTARl ,  XSTARl 

[118]  r  AND  ESTAR1 

[119]  r TRANSFORM Y 

[120]  YSTARn<rPn+  .xYSTARl 

[121]  o  TRANSFORM  X 

[122]  XSTARH+PU+ .xXSTARl 

[123]  r  TRANSFORM  RESIDUALS 

[124]  ESTAR^<rPH  +  .  xESTARl 

[12  5]  r  CALCULATE  VALUE  FOR  2- STEP  BETA 

[126]  S2S-e-yS2’^fl4@XS2,.4i?4 

[127]  B2SIZC0UNTI1+B2SZ1\1 

[128]  B2SS[COfWZ7J]«-B2S[2;] 

[129]  r  CALCULATE  THE  SCHMIDT  STATISTIC  BEFORE 

[130]  r TRANSFORMATION 

[131]  DS1C+ ( (§EH AT )+. xEHAT) 

[132]  DS+,(.(§EHAT)  +  .x(AS+.xEHAT))*DS1C 

[13  3]  r  COMPARE  TABLED  VALUE  FOR  SCHMIDT  STATISTIC 
[134]  r  TO  ABOVE  CALCULATED  VALUE  FOR  SCHMIDT  STATISTIC 
[13  5]  DEBSMT<-(DS£DSMTL  ) 

[136]  o  IF  DETECTED  INCREMENT  1 

[137]  DSBCNT+DSBCNT+DEBSMT 
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[138]  INSBCT+ ( ( DSMTL<DS ) a ( DSMTU>DS ) ) 

[13  9]  fi  IF  INCONCLUSIVE  INCREMENT  1 

[140]  STIBCT+STIBCT+INSBCT 

[141]  o  USE  1,4  P  MATRIX  TO  TRANSFORM  Y ,  X ,  AND  EH  AT 

[142]  fi  TRANSFORM Y 

[143]  YSr4tfl4-*-P14  +  .xy 

[144]  fi  TRANSFORM  X 

[145]  XSTARm+Pm  +  .xX 

[146]  fi  CALCULATE  VALUE  FOR  1,4  BETA 

[147]  Bl4«-YS771P14gX,ST4P14 

[148]  BimZCOUNTIl+Bmil;l 

[149]  BmSZCOUNTI]*-Bmi2 ;  ] 

[150]  r  CALCULATE  RESIDUALS 

[151]  ESTARm+YSTARlU- (XSTARin+.  xfll4  ) 

[15  2]  fi  CALCULATE  THE  SCHMIDT  STATISTIC  AFTER 
[15  3]  fi  THE  TRANSFORMATION 

[154]  DS2D+((§ESTARin)+ .xESTARlV) 

[155]  DS+,  ((§ESTARm)+.x(AS+.*ESTARin))iDS2D 

[156]  fi  COMPARE  TABLED  VALUE  FOR  SCHMIDT  STATISTIC 

[15  7]  fi  TO  ABOVE  CALCULATED  VALUE  FOR  SCHMIDT  STATISTIC 

[158]  DETSMA+(DS<DSMTL) 

[159]  fi  IF  DETECTED  INCREMENT  1 

[160]  DSACNT+DSACNT+DETSMA 

[161]  INS  ACT*-  (  (DSMTL<DS  )a  (DSMTU>DS  )  ) 

[16  2]  fi  IF  INCONCLUSIVE  INCREMENT  1 
[163]  STANCT+STANCT+INSACT 

[16  4]  ->  (COUNTI<5  ) /INNER 

[16  5]  fi  CALCULATE  AVERAGE  VALUES  FOR  BETAS 

[166]  BLSIAVlC0UNT01+(+/BLSI)*5 

[167]  BLSSAVZC0UNT01+(+/BLSS)* 5 

[168]  B2SIAVZCOUNTO]*-(+/B2SI)*5 

[169]  B2SSAVZC0UNT01+(+/B2SS)*5 

[170]  BmiAVZC0UNT01*-(+/BlUI)i5 

[171]  BmSAVZCOUNTO]  +  (+/BmS)*5 

[172]  +(COUNTO<VO)/ OUTER 


73 


[173]  AF+DSACNT 

[174]  BE+DSBCNT 
[17  5]  LSI+BLSIAV 

[176]  LSS<-BLSSAV 

[177]  I2S<rB2SIAV 

[178]  S2S<-F2SSA7 

[179]  IlU+BinlAV 

[180]  S14«-B14SA7 
V 
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